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1 . INTRODUCTION 
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1 . 1 Summary 

This report details research performed under NASA Grant NAG3-491 , 
with startup date of November 15, 1983 and completion on December 31, 
1985. The research centers on evaluation of thermoviscoplastic 
constitutive models for metals at elevated temperatures. The primary 
intent of the research was threefold in nature: 1) to improve existing 
experimental facilities within the Mechanics and Materials Labs at Texas 
A&M University in order to perform the complex experiments required 
under the grant; 2) to compare existing models to experiment for the 
materials to be discussed herein; and 3) to extend existing models where 
necessary to better predict complex material response. The research 
group consists of Drs. W.E. Haisler and D.H. Allen, assisted by five 
graduate research assistants. 


1.2 Acknowledgement 

The authors express their thanks for the support provided for this 
research by the NASA Lewis Research Center. The technical grant monitor 
is Dr. R.L. Thompson. 


2. RESEARCH OUTLINE 
2.1 Summary of Completed Research 

The important details of the research are outlined in this 
section. The findings are divided into the following six sections: 1) 
numerical integration techniques; 2) thermodynamics and internal state 
variables; 3) experimental lab development; H) comparison of models at 
room temperature; 5) comparison of models at elevated temperature; and 
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o) integrated software development. These topics are summarized in the 
following seven sections (2.2 through 2.7). Further details are given 
in the technical reports in the appendix. A summary of important 


findings is given in Section 2.8. 


N86- 30228 X>, - 3; 


2.2 Numerical Integration Techniques J7 7 

The sensitivity of the unified constitutive theories proposed by 
Bodner, Walker, Krieg, and Miller, to numerical integration techniques 
and slight changes in material parameters was investigated. Evaluations 
were based upon numerical simulations of Hastelloy-x at 1800°F, in 
uniaxial form, by specifying input strain histories and comparing output 
stress histories. The constitutive models selected for this study were 
chosen primarily because of the availability of a prescribed methodology 
for material parameter evaluation. In addition, they appear to be the 
most qualitatively attractive theories available and thus warrant 
further evaluation. 


The investigation begins by presenting the various models in 
uniaxial differential equation form, followed by a description of the 
methods used for material parameter evaluation. The main thrust of the 
research was to assess the degree of sensitivity of each model to 
integration techniques and/or variations in the respective material 
constants. An example is shown in Fig. 1 for Bodner's model. 

For the numerical time integration study, four commonly used 
algorithms were selected. These included Explicit Euler Forward 
Difference, Implicit Trapezoidal, 4th Order Runge-Kutta, and Trapezoidal 
Predictoi — Correction methods. In addition, four different strain rate 
input histories were selected to insure that evaluations were not biased 
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Fig. 1 . Comparison of Integration Methods For Bodner’s Model 
With Equivalent Computation Time Allowed For Each 
Method (Hastelloy-X at 1800° F). 


to one specific simulation. The strain histories considered were not 
developed to model any specific physical phenomena. However, they are 
representative of conditions under which a material could be strained 
either in the laboratory or the f’ield. 

The integration investigation was divided into two tasks. The 
first task included a comparison of the numerical stability of each 
constitutive theory in order to assess the degree of nonlinearity of the 
prescribed growth laws. The comparisons were based upon predicted 
stress-strain behavior for constant strain rate tests simulating 
monotonic loading. Explicit Euler Forward integration was used to 
obtain stress histories for various time steps. It was concluded that 
Miller's model was the most sensitive to step size variation because of 
the oscillatory response of the stress-strain curve for the larger At 
values. It is postulated that the oscillatory nature of Miller's model 
can be traced back to the hyperbolic sine function that is used to 
characterize the inelastic strain rate and other growth laws. 

In addition, it was found that the predicted stress response for 
each model appears to be "self correcting" for large time increments 
when the constitutive equations are integrated by the Euler method. 
This "self-correcting" phenomenon appears to be an intrinsic property of 
the prescribed growth laws and the coupled nature of the equations. 

The second task of the integration investigation was to evaluate 
the constitutive theories in terms of solution stability, accuracy, and 
computational efficiency when numerically integrated using the various 
algorithms considered. A test matrix consisting of integration method, 
time step sizes, and strain input histories was established in order to 
make qualitative and quantitative comparisons. For each strain history 
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considered, a baseline or "pseudo correct" solution was obtained 


analytically, since no experimental data were available. 

The computational efficiency of each integration method was 
obtained through selection of equivalent At’s. Euler’s method was 
observed to be the most efficient, followed by the Implicit Trapezoidal 
method, then the 4th order Runge-Kutta method, and finally the 
Trapezoidal Predictoi — Corrector method. In addition, each model was 
evaluated for its computational efficiency. The results indicate that 
Walker's model resulted in the fastest execution time, followed by 
Krieg, Bodner, and Miller, respectively. 

In terms of qualitative and quantitative comparisons the following 
observations were made: 1) In general, for equal time step size and 
equivalent computation integration times, the 4th order Runge-Kutta 
method consistently overpredicts the state of stress while the 
Trapezoidal method generally underestimates the stress history. The 
Euler and Predictoi — corrector methods appear to provide accurate results 
at any point in the simulation; 2) the amount of stress-over3hoot (for 
the larger step sizes) appears to grow with increasing strain rate; and 
3) no "self-correct ing" tendency was observed when either the 4th order 
Runge-Kutta or Trapezoidal method was used. 

Another important issue this research addressed was what amount of 
degradation in predicted response could be expected due to variations in 
material constants. Two procedures were developed to study this. The 
first procedure characterized the sensitivity of each constant within a 
particular constitutive model. Each constant's sensitivity to variation 
was determined by adjusting it by 5 % and comparing predicted results. 
The second procedure provided an "upper bound" on the effect that 
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experimental error would have on the determination of the material 


parameters and ultimately, predicted response. The effect of 
experimental uncertainty was studied by adjusting numerical test data by 
±5? and then using the adjusted data to determine new material 
constants. 

From the results of the sensitivity study, the following 
observations were made: 1) For the load history considered, 
constitutive models can have up to a 12$ variation in predicted response 
for as little as a 5$ variation in a given material constant; and 2) 
Miller's model was shown to be the least sensitive to material parameter 
variation possibly due to its strong microphysical basis. On the other 
hand, Walker's model was the most sensitive, showing a 30$ over- 
prediction in stress during the relaxation period of the load history 
considered. Bodner's and Krieg's models were found to be between these 
two extremes. Further details on this phase of the research can be 
found in Appendices 6.1 and 6.2. 

It should be noted that the conclusions reached from this study 
are based on the use of uniaxial simulations. It is expected that the 
conclusions would hold for multiaxial situations. However, further 
evaluation is needed in this area. 

2.3 Thermodynamics and Internal State Variables 

As a part of the research effort it was concluded that a casting 
of the currently available constitutive models into a common framework 
would help in identifying the relative merits and/or deficiencies of 
each model. For this purpose it was decided to cast the models into the 
framework of internal state variable (ISV) theory, in which a set of 
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state variables which are not observable are utilized to account for 
thermodynamically dissipative processe s s uch as dislocation movement, 
grain boundary sliding, deformation twining, phase changes, etc. This 
framework has been used to determine qualitative differences in the 
models studied under the current grant. Furthermore, in this process it 
was found that certain thermodynamic issues could be clarified. 

The specific ISV framework uses the approach initially proposed by 
Coleman and Gurtin [1], in which the observable state variables are 
supplemented with a set of internal state variables: 

a ij = a ij ( V t)( K = 1>2> » U> 


where are a set of n internal state variables, in this case 
designated as second order tensor functions of space and time t. 

The method of Coleman and Noll [2] may be utilized to obtain 
thermodynamic constraints on the class above of materials including a 
broad range of thermoviscoplastic metals. The results are: 
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where o, , is the stress tensor, e, . is the infinitesimal strain tensor, 
kl kl 

p is the mass density, and h is the specific Helmnoltz free energy, 
given by 


h = u-Ts = h( e , 
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where u is the internal energy, T is the temperature, and s is the 
entropy. The Helmholtz free energy is then expanded in a Taylor series 
in strain, temperature, and the ISV's and substituted into (2) to obtain 
the stress-strain relation: 


*ij °ijkl (e kl 


- e 


kl 


“kl 


(4) 


T 

where C. .. . is the elastic modulus tensor, and e , is the thermal strain 
ljkl kl 

tensor. The above equations serve as a framework of comparison for all 
of the models considered in this research. These are adjoined with ISV 
growth laws of the general form: 


k n k f 

a. .= fl . . ( e , 
lj xj mn’ 
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a ) 
mn 


(5) 


The ISV's represent locally averaged measures of dislocation arrangement 
(back stress) dislocation density (drag stress), etc. ISV growth laws 
(5) represent the principal differences between the models studied in 
the current research effort. These equations have been itemized in 
Table 1 for ten different current models. Qualitative differences in 
the models are also discussed in the table. 

The main content of this part of the research has been to review 
and clarify the continuum and thermodynamics based internal state 
variable model for application to thermoviscoplastic metals. In this 
process the following points have been made: 

1) the definition of an internal state variable utilized in this 


model has been clarified; 
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2) internal state variables in metals represent local averages of 
dislocation arrangement, dislocation density, and intergranular damage, 

3) in the context of the ISV definition given here, inelastic 
strain may also be interpreted as an internal state variable; 

M) the path dependent Helmholtz free energy may be expanded in a 
second order expansion in elastic strain and temperature in order to 
obtain a stress-strain equation of state; 

5) rate dependence enters the constitutive equations implicitly 
via the inelastic strain, as demonstrated by the nonlinear standard 
solid analog; and 

6) a three-dimensional generalization of the standard solid may 
be used as a means of comparison of the general form of several 
currently proposed models. 

Further details of this development are contained in Appendix 6.3. 
2. A Experimental Lab Development 

This section describes the test equipment used to carry out the 
tests required for this work. The first sub-section describes the test 
set-up used for the creep and strain transient tests based on a dead 
weight creep frame. Topics to be covered include the load frame, data 
acquisition, and temperature* measurement and control. The next sub- 
section describes the test set-up used for the monotonic strain rate, 
fully reversed cyclic tests, and transient tests during cyclic loading 
as well as creep and strain transient tests based on a computer 
controlled test set-up. 
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Creep Frame Test Set-up 


Load Frame , The load frame utilized was a creep frame produced at 
Texas A&M University. The initial configuration included a constant 
load cam and a constant stress cam as described by Garofalo, et al. [ 3 ]. 
The lift for the pan was provided by an ATS (Applied Test Systems) 2081 
cyclic load module. Several modifications were made to this set-up. The 
constant stress cam was removed and a second constant load cam 
installed. This increased the maximum weight limit from 400 lbs. to 800 
lbs. This also increased the stability of the loadpan during load-up 
and provided a seven to one load magnification. The load pan was 
increased in size and supports were added to further aid stability 
during load-up. The cyclic load table was removed and replaced with a 
Century-Fox model CF-59 5 ton capacity hydraulic jack. This provided 
more lift capacity and gave more room for machine deflection during 
load-up. The p inal configuration of the machine is shown in Fig. 2. 

Further improvements included a more advanced load pan which would 
provide automated removal and replacement of portions of the load. 
Jones, et al. [4] have described such an apparatus. This modification 
is needed to improve the results of the strain transient tests. A 
second improvement was a more advanced method of load-up. For this 
purpose, a pneumatic jack was included. 

Data Acquisition . The load was measured with a Strain sert TLN20- 
256K Tension Link driven by a Hewlett Packard 8805A Carrier Preamp. The 
strain was measured with an ATS model 4112 LVDT (Linear Variable 
Displacement Transducer) and the extensometer was driven by an ATS model 
6974 signal conditioner. Two Hewlett Packard 8803A Low Level Preamps 
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Fig. 2. Creep Frame Final Configuration 








and two low pass filters were necessary to achieve the necessary 
resolution. Load and strain measurements were recorded on a Hewlett 
Packard series 7700 strip chart recorder. 

This data acquisition system was very prone to drift, noise, and 
nonlinearity. All strain amplification, filtration, and data 
acquisition equipment were eventually replaced with a Kiethly model 197 
logging digital multimeter. This set-up proved much more stable. 
However, it lacked the ability to take in a large number of data points. 
Much improvement is needed in strain data acquisition to produce a 
system which is stable, can take in large numbers of points, and has 
sufficient resolution. A common set of transducers and a common data 
acquisition system should be utilized for both load frames. This would 
negate errors arising from using two different systems. 

Temperature Measurement and Control . An ATS 2961 Clamshell oven 
and an LFE series 230 temperature controller were used for temperature 
control. Temperature measurement was handled by two 20 gauge K-type 
thermocouples. These were placed in contact with the sample in the 
middle of the gauge section on diametrically opposite sides. One 
thermocouple was used as input to the temperature controller. The other 
thermocouple was input to a Fluke 2176A digital thermometer for readout. 

Improvements in this system would include the ability to monitor 
the temperature at several points simultaneously along the gauge 
section. The temperature should also be input to a data acquisition 
system. The method of thermocouple attachment should also be upgraded. 
The optimum method of thermocouple attachment would be individual 
welding of the leads to the sample surface. This is described in the 
ASTM Thermocouple Handbook [5]. Such a procedure was attempted with 20 
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gauge thermocouple wire and a Duracom thermocouple welder. The samples 

tended to fail prematurely at the . thermocouple welds. Therefore, this 

method of thermocouple attachment was abandoned. Another possible 
method using smaller 28 gauge thermocouples was used with the strain- 
controlled set-up and will be discussed in that section. A common 
temperature measurement and control system for both load frames would be 
another improvement. This would eliminate relative errors between the 
two systems. 

Computer Controlled Test Set-up 

Load Frame . The load frame utilized in this set-up was an MTS 
(Materials Test System) model 880 electro hydraulic testing machine. 
MTS 652.01, water-cooled, hydraulic grips allowed fully : reversed cyclic 
tests to be carried out at high temperature. The frame was controlled 
by a Digital Micro PDP-11. Computer programs were written to run 
monotonic tension tests, cyclic tests, cyclic tests with hold times, 
creep tests, and creep stress drop tests. The load frame is shown in 
Fig. 3. 

Data Acquisition . The Micro PDP-11 also handled data acquisition 
functions. An MTS 661.21A-02 10 kip load cell was the load transducer. 
An MTS 632.41B-02 axial extensometer was the strain transducer. This 
device had quartz extension rods which contacted the sample at two 120° 
punch holes. All data were stored on 5.25 inch floppy diskettes and 
could be retrieved in hard copy or graphical form. 

Temperature Measurement and Control . An MTS 652, three-zone, 
clamshell furnace and three Research Incorporated 63911 Process 
Temperature and Power Controllers were used for temperature control. 
Temperature Measurement was handled by six 28 gauge K-type 
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thermocouples. These were placed in contact with the sample. Three 
thermocouples were fed into a Fluke 2176 a Digital Thermometer for 
readout. These were placed with one each at the top, middle, and bottom 
of the gauge section. The other three thermocouples were fed into the 
temperature controllers. These were placed in the center of the furnace 
zone each was to sense. One thermocouple was placed in the center of 
the gauge section and one on each grip. 

The thermocouples were supported at the grips by fiberglass 
thread. They were fixed to the sample by self-supporting means. The 
thermocouples at the top and bottom of the gauge section were wound 
around the sample. The thermocouples used in the center of the gauge 
section were brought into the oven from different directions and tied to 
each other. These thermocouples were then wound around the sample. The 
thermocouples used in the center of the gauge section were brought into 
the oven from different directions and tied to each other. These 
thermocouples were then wound around the sample for contact. Welding 
the thermocouples to the sample would have produced harder contacts with 
more reliable temperature measurement. However, as mentioned earlier, 
premature failure occurred at the welds. 

Further details about the experimental development are given in 
Appendix 6.5. 

2.5 Comparison of Models at Room Temperature 

The purpose of the room temperature testing was to compare three 
theories based, on considerations of microstructural behavior. These 
theories are those of Krieg, et al., Bodner, et al., and Miller, et al. 
(references [18], [107], and [115], respectively, in Appendix 6.4). 
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Each theory is proposed in uniaxial form to simplify experimental and 
analytical analysis. Since each theory is proposed for a different 
material, a set of experiments is performed with a candidate material in 
order to create a common experimental data base with which the 
mathematical representation of each model can be characterized. These 
theories are subsequently analyzed qualitatively and quantitatively over 
a broad range of material behavior to determine their suitability for 
predicting inelastic deformation. 

The experimental data base requirements for each theory are as 
follows. The theory of Bodner, et al., requires several constant strain 
rate tensile tests and several creep tests. The theory of Miller, 
et.al., requires several creep tests, several constant strain rate 
tensile tests, and some information about the material such as the 
melting point and the activation energy. The theory of Krieg, et al., 
requires a constant strain rate tensile test, a creep curve complete 
with primary and secondary regions, and several stress-drop tests in 
which the stress is reduced rapidly during steady-state creep in order 
to examine the resulting transient strain response. In each case the 
actual number of tests, unless specified, is dependent on the amount of 
data deemed necessary for accurate evaluation of the material constants. 
The candidate material chosen is the aluminum alloy 5086 in the form of 
uniaxial bars, and the material tests are performed at room temperature. 

In determining the material constants for each theory from these 
tests, the procedure proposed by the authors is followed. It is found 
that some of these procedures may be difficult to implement. In 
particular, difficulties are encountered in interpreting the stress-drop 
test, in which the quantity of interest is a zero strain rate or merely 
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the absence of a resolvable strain rate a^ter unloading. It can be seen 
that the results of this test are somewhat dependent* on the resolution 
of the experimental equipment. In addition, the values of several 
constants in the theories are set arbitrarily, possibly leading to 
inaccuracies in the ensuing predictions. Many procedures require 
extensive data reduction and graphing, making an interactive data 
acquisition and computer system an indispensable tool. Finally, 
correction of the value of any one constant can be a difficult task, as 
some of the constants of each theory may be interdependent. 

Each constitutive theory is presented by highly nonlinear, 
numerically "stiff" differential equations. Since the objective of this 
research is not to determine numerically efficient integration 
techniques for these equations, a stable first order forward integration 
scheme was used. 

Numerical simulations were performed to compare the predictive 
capability of each theory to experimental results for constant strain 
rate tensile tests, constant strain rate cyclic tests, and complex 
strain rate history tests. As an example, Fig. 4 compares theory to 
experiment for the first cycle of a constant strain rate cyclic test. 
The experiment is modelled fairly well, but it can be seen that the 
elastic-inelastic transitions are too sharp, possibly indicating a need 
for more complex hardening laws. In addition, the Bauschinger effect 
may not be modelled well by Bodner's theory and Miller's theory due to 
the lack of a representation for the back stress. The discrepancy in 
the width and height of the hysteresis loop predicted by Krieg's theory 
may be due to the ambiguous measurements of the back stress in the 

stress-drop tests, leading to possible inaccuracies in the material 
constants. Results and conclusions of comparisons of other tests may be 

found in detail in Appendix 6.5 of this report. 
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Fig. 4. Comparison of Models to Experiment for IN 718 at 1100 F 
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2.6 Comparison of Models at Elevated Temperature 

The theories used in this work Include Bodner's anisotropic model [6]; 
Krieg, Swearengen, and Rohde's model [7]; Schmidt and Miller's model [8]; and 
Walker's exponential model [6]. These models have been chosen because they 
are widely publicized, examples of constant calculation methods have been 
presented, and the attempt is being made to extend them to non- isothermal 
modelling. 

Bodner's anisotropic model utilizes an inelastic strain rate equation 
extended from classical plasticity theory. The inelastic strain is assumed to 
be exponential in stress. A directional parameter has been included as an 
incrementally isotropic variable. The growth laws are microphenomenological ly 
based and model hardening, dynamic recovery, and static thermal recovery. 
Bodner's model is given below: 


.D = ^ D o exp [- \ (^) 2n ] sgno 

Z = Z l + Z A = Z 1 + B sgno l 
i l = m l( Z r Z 1 ) W p - A x Z x (— 

B = M 2 (Z 3 sgno-Z A ) W - A 2 Z^-^) 


sgn Z A 


( 6 ) 


Krieg, Swearengen, and Rhode's model is based on microphenomenological 
considerations. The inelastic strain rate equation is based on a power law 
formulation. The growth laws model hardening and static thermal recovery. 
This model is given below: 


E 1 - C <^) n sgn (o-B) 

I 2 A 3< B - B o )2 

8 - Aj s - A 2 { 8-B 0 ) (e -1) sgn(B) 
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(7) 


0 - *41**1 - A 5 (°- |) 0 )n 


Schmidt and Miller's model is a microphenomenologically based model which 
uses a hyperbolic sine inelastic strain rate equation. The growth laws for 
back stress and drag stress model hardening and static thermal recovery. 
Parameters have been added to this model which account for strengthening due 
to solutes. Schmidt and Miller's model is given below: 


i 1 - B ‘ sinh[(-|r“ s 9 n (f -B) 

sol 

B = H, i 1 - HjB ' sinh (A^B)) n sgn(B) 

D = H 2 ! e 1 1 (C 2 + |B | - ^ D 3 ) -H 2 C 2 b' sinh (A^ 1 * 5 ) n 


F sol F sol max 


exp (. !l°.g.(. lp .l..k.io g .(J) )) 2 


( 8 ) 


Walker's model has an exponentially based inelastic strain rate 
equation. The back stress growth law models hardening dynamic recovery, and 
static thermal recovery. The dynamic recovery term has been modified to 
handle negative strain rate sensitivity. The drag stress variable allows for 
cyclic hardening or softening. Walker's model is presented below: 

.1 . exp (°-B) - l sgn(o . B) 

B = n 2 ^-B [n 3 + n 4 exp ( -n 5 1 1 og (— ) 1 1 R + n g 

*o 

D = D^+ D 2 exp (-n^R) 

R = II 1 ! (9) 


- 21 - 



2.6.1 Material Response 

The material used in this work is Inconel 718 at 1100° F. All samples 
have been sent through a common heat treatment prior to being tested. This 
material system has several distinctive characteristics. The material 
responds with negative strain rate sensitivity. Discontinuous yielding or the 
Portevin-Le Chatelier effect is observed. The material cyclically work- 
softens and has a fatigue life of 5 to 30 cycles when cycled at strain limits 
over ± 1 percent strain. The yield point of the material is in the region of 
95 to 110 ksi. 

2.6.2 Experimental Considerations 

The material samples were configured to meet ASTM standard E606-77T. The 
tests were carried out on a computer-controlled, MTS 880, 110 kip, electro- 
hydraulic testing machine. A Digital Micro PDP-11 computer controlled the 
tests and acquired the data. An MTS 652 furnace system with three-zone 
clamshell furnace and hydraulic grips was utilized. An MTS 632-41-02 high- 
temperature axial extensometer was utilized for displacement measurement. 
Temperature was monitored with K-type thermocouples and a Fluke 2176A digital 
thermometer. Welding of the thermocouples to the samples initiated failure 
and a self-supporting, surface contact method was used for thermocouple 
attachment. 

The test program consisted of the following tests: 

(a) 2 monotonic tension tests to 1.556 strain 
([1] 3.15 xl0 -3 sec _1 [2] and 7.25xlO" 6 sec” 1 ) 

(b) 5 fully reversed cyclic tests to ± .8% strain 

-3 -1 -6 -1 

(strain rates between 1x10 sec and 7x10 sec” ) 
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(c) 5 constant load creep tests 

(applied stresses between 120 ksi & 140 ksi) 

(d) Back stress measuring tests during cyclic loading and during secondary 
creep 

(e) a complex history test 
2.6.3 Determination of Parameters 

Beek [4] has shown that the parameter calculation process presented with 
each model can produce constants which may produce undesirable model 
response. A method for determining material parameters which lessened this 
problem and provided some commonality in the calculations was used. A set of 
initial assumptions was made based on the expected response of the material. 
These initial assumptions included the following: 

(1) Back stress was responsible for hardening in monotonic tension; 

(2) Drag stress was responsible for cyclic softening; 

(3) Thermal recovery could be neglected for rapid tests 
(e T > lxl0 -4 sec -1 ) 

(4) Drag stress thermal recovery could be calculated from slow cyclic tests; 
and 

(5) Back stress thermal recovery could be calculated from creep tests. 

Hand calculations based on these assumptions provided initial estimates 
of the material parameters. Computer-aided iterations were then used to tune 
the response to match specific data points. The parameters which resulted 
from the iteration process indicated that the following initial assumptions 
would have been more appropriate: 

(1) the inelastic strain rate equation constant could be set to provide the 
proper scaling and strain rate sensitivity; 
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(2) back stress hardening produced monotonic hardening; 

(3) solution strengthening parameters produced negative strain rate 
sensitivity; 

(4) drag stress hardening could be used to set cyclic work softening 
characteristics; and 

(5) thermal recovery effects were small and masked by solute strengthening. 

A comparison of the response of the solute strengthening corrected models 
to experiment is shown in Fig. 5. These results provide several conclusions 
concerning various parts of the models. The theories of Walker and Bodner 
with exponentially based inelastic strain rate equations handle the negative 
strain rate sensitivity most effectively. The theories of Krieg, et al.. 
Miller, and Bodner produce oversquare stress-strain curves. 
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The addition of dynamic recovery terms to the directional growth law 
appears to help the strain rate sensitivity problem. Bodner's corrected model 
shows less response to strain rate jumps. This is due to hardening based on 
plastic work rate. Walker's correction for solute strengthening comes closest 
to reproducing the strain rate sensitivity. The correction applied to 
Bodner's model produces unstable response and a correction similar to Walker's 
should be considered. The correction factors negated the effects of thermal 
recovery during hold times for Bodner's and Krieg, et al's. models. The drag 
stress growth law of Walker provided the closet fit to data over several 
cycles at the higher strain rates. Lower strain rates showed Bodner's model 
handling cyclic response most effectively. Walker's model suffers from the 
lack of drag stress thermal recovery at the lowest strain rate. 

These results show that the combi nation of exponential strain rate 
equations with dynamic recovery work best for this material system. Isotropic 
thermal recovery and exponential growth of the isotropic internal state 
variable are useful for modelling cyclic response. 

Further details of this portion of the research are given in Appendix 

6.5. 


2.7 Integrated Software Development 

Because of the vast number of constitutive theories available, only four 
candidate models were selected for this research. These include the work of 
Walker[6,9], Krieg, Swearengen, and Rohde[7], Bodner[6], and Schmidt and 
Mi Her [101. These models were chosen primarily because they appear to be the 
most qualitatively attractive theories available and thus warrant further 
evaluation. For this paper, only the procedure developed for Walker's model 
will be presented. 
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The model proposed by Walker is a viscoplastic theory which uses an 
exponential type inelastic strain rate relation. The growth law modeling back 
stress is of the hardening/recovery format and accounts for both dynamic and 
static thermal recovery. The drag stress term models isotropic hardening, 
thus taking into account the cyclic hardening or softening characteristics of 
a material. The uniaxial differential form of Walker's exponential model may 
be written as: 


exp(- 


-B 


e = 


_D_ 

6 


) - 1 


sgn(o - B) 


B 

D 


• I R 

n 2 e - B(ln 3 + n„ exp(-n 5 j log(— ) ] R + n 6 > , 

Dj - D 2 exp(-n 7 R) , R ° 


( 10 ) 

( 11 ) 

( 12 ) 


R - |e'| , ( 13 ) 

where a is the applied stress, e 1 is the inelastic strain, B is the back 
stress, and D is the drag stress. A superposed dot above the variables 
denotes differentiation with respect to time. In addition, 6 , n2» 03, n^, 
n 5 ’ n 6» n 7» °2 and °2 are ma terial parameters. Therefore, this model requires 
nine constants to be evaluated, along with selecting Young's Modulus E and the 
strain aging parameter . 

The tests required to determine the constants for Walker's model using 
the following procedure include: 1 ) A series of constant strain rate steady 
state hysteresis loops under fully reversed strain controlled conditions; 2) 
cyclic hold tests performed on the unloading branch of the cyclic tests; and 
3 ) long term monotonic tension tests. The monotonic tension tests may not be 
necessary if acceptable values of the limiting stress 0^ can be obtained 
from the first half cycle of the cyclic tests. 
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versus 


The determination of the constants begins by plotting o^ m 

In (e*) . A nonlinear representation signifies that strain aging and/or 

thermal recovery effects are present and thus need to be modeled. If is 

not obtained experimentally, it can be estimated in a manner similar to that 
proposed by Lindholm, et al.[ 6 !. 

Under conditions of uniaxial tension loading, when e 1 is assumed to be a 
constant and equal to the applied strain rate, eq. ( 10 ) may be written as: 

o = D In (Be 1 ) + B . (14) 

Using the evolution equation defining the back stress (eq. (11)) and on the 
physical basis that D remains constant during monotonic loading, 

then do/de 1 (or o) may be written as: 

o = n 2 - bI n 3 + i\exp(-n 5 1 log (— 7 —) | ) + r^/e 1 ] . (15) 

Thus, equations (14) and (15), can be combined yielding 

O = -No + [n 2 + NDln(Be 1 )] , (16) 

where 

N = n 3 + n 4 exp (-n 5 1 log (- 7 ^) | ) + ng/e 1 . (17) 

R 0 

Therefore, equation (16) indicates that a plot of 0 versus o will be linear, 
having a slope of N and an x-intercept of o^ m . 

The constant n$ is computed by determining where the effect of strain 
aging is considered negligible and may be written as: 
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( 18 ) 


n 5 = - ( 1 n( t ) / 1 log (-?-) | ) . 

R o 

The constants R 0 and represent the strain rate at which the strain aging 
correction is a maximum and minimum, respectively, and x denotes the residual 
correction at rate R : . It should be noted that x also affects the rate of 

decay of the strain aging correction and selection of too small of a value 
will result in a very localized correction. 

The next step in this procedure is to compute the dynamic and static 
thermal recovery constants ng, n^, and ng using equation (17). Equation (17) 
can be rewritten a number of times, corresponding to the different monotonic 
tests as: 

N i = n 3 + n l) f i + n 6 /e| , (19) 

where 

f. = exp(-n 5 |log (~~-)| ) . (20) 

R o 

Thus, the three parameters, ng, n 4 , and ng, may be obtained simultaneously 
using a multiple linear regression scheme. 

The constant n 2 is computed on the basis that B saturates to B^ m at 
large inelastic strains. Hence, B=0 and equation (11) reduces to 


lim 


n, + n f + n,/£ 


•I 


( 21 ) 


If it is assumed that the ratio /B„ „ will remain the same for the 

exp exp 

limiting condition at sufficiently large inelastic strains, then 


°exp 

B exp 


q lim 

B lim 


( 22 ) 
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Substituting eq. (21) into (22) and solving for n£ results in 
In, + n„f ♦ n 6 /c>]„,. m B 

n = C_ 


( 23 ) 


exp 


where a-^ values are obtained from long term monotonic tension tests 
or e-plots and o gxp and B gxp values come from cyclic hold tests. 

The initial value drag stress Dq and the inelastic strain rate scaler b 
are determined by rewriting equation (14) using the limiting values of o and 
B as: 


° 11 m - B lim ■ O o'n<;’> + Ooln(s) • (Z«) 

Since B lim is given by equation (21), it can be substituted into (24) 
resulting in 


6 oB = D o 1n (e ) + D o 1n ( e ) 

where 


6 oB ' °lim 


+ nj 


•I 


(25) 

(26) 


Equation (25) indicates that a plot of 6 aB versus ^(e 1 ) should be linear, 
having a slope of Oq and an intercept of D ln(e) from which b can be found 
directly. 

Up to this point, the only tests that were needed, in order to compute 
the material parameters, have been monotonic tension and cyclic hold tests. 
To obtain the isotropic hardening and recovery constants D^, D 2 . and n^, 
saturated cyclic hysteresis data are required. By estimating the cumulative 
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inelastic strain from applied stress, E, and strain amplitude and assuming 
that D saturates to Dj, then n-/ can be approximated as: 


n 7 = -ln(t)/R avg , 


(27) 


where R avg is the average of R for a number of tests and t is an approximation 
of zero. 

On the physical basis that B saturates much more rapidly than D, equation 
(14) can be written as: 




where 


°lim B 1 i m 

1 n (ge 1 ) 


B lim 


n 2 /N . 


(28) 

(29) 


The values of o lim and N in equations (28) and (29) are estimated from the 
o-plot after cyclic saturation has occurred. If several tests are used, 
would be computed as the arithmetic mean. 

The final parameter to determine is D 2 . When R=0 equation (12) reduces 
to 


D 0 = D 1 - D 2 , (30) 

from which D 2 may be computed directly, completing the constant calculation 
procedure. 

The computer algorithm used to compute the material parameter proceeds in 
the following manner: 

1) Values of o^- m and N from monotonic tension tests are computed using 
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equation (16) and a least squares procedure. 

2) After selecting R 0 , R, and t, n 5 is evaluated using equation (18). 

3) The parameters n^, n^, ng are determined from equations (19) and (20) 
using a multiple linear regression scheme. 

4) The back stress hardening coefficient n 2 is computed directly by 
equation (23). 

5) The constants Dq and e are evaluated using equations (25) and (26) 
via a least squares procedure. 

6) After computing the cumulative inelastic strain from the experimental 
data, n 7 is evaluated directly from equation (27). 

7) Dj and D 2 are computed using equation (28) and (30). 

A method for obtaining the material parameters for Walker's model has been 
presented which is a synthesis of both physical and numerical 
approximations. The associated computer algorithm allows the user to use 
either a totally automated procedure or engineering intuition at selected 
points when computing constants. Numerical simulations of Inconel 718 at 
1100°F are to be presented. 


2.8 Conclusion 

Several general conclusions can be made from the current research. These 
are as follows: 

1) Euler's method is the most efficient for integrating the models of 
those methods considered herein; 

2) Walker's model is the most economical to integrate, while Miller's 
model is least economical and also tends to be unstable; 

3) 1SV theory can be used as general framework for comparing all of the 
models in the current literataure; 

4) a laboratory has been developed at Texas A&M which is capable of 
performing all of the complicated elevated temperature experiments necessary 
to perform the research detailed herein; 

5) for aluminum at room temperature, none of the models considered 
herein appears to be very accurate in predicting uniaxial cyclic response; 

6) at 1100°F Walker's model appears to be most accurate for predicting 
the uniaxial cyclic response of IN 718. 
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Several unified constitutive models were tested in uniaxial form by 
specifying input strain histories and comparing output stress histories. 
The purpose of the tests was to evaluate several time integration methods 
with regard to accuracy, stability, and computational economy. The sen- 
sitivity of the models to slight changes in input constants was also in- 
vestigated. Results are presented for IN1Q0 at 1350°F and Hastelloy-X at 
1800°F. 


INTRODUCTION 


The characterization of the constitutive behaviour of metals has its 
roots in the early work of Tresca, Levy, vonMises, Hencky, Prandtl, Reuss, 

Prager, and Ziegler (Refs. 1-8). These early models are incremental in 
‘ nature, assume that plasticity and creep can be separated, and they incor- 
porate a yield function, flow rule, and hardening rule to define the plastic 
strain increment. These original incremental theories have been expanded 
and modified by many researchers so that they provide adequate, and often 
very good predictions of rate- independent plastic flow (see for example 
Refs. 9-10). However, they are sometimes criticized as having no formal 
micromechanical basis upon which to make the assumption of an uncoupling %, 

of the inelastic strain into rate-independent (plastic) and rate-dependent 
(creep) strain components. Nevertheless, the classical incremental theories 
are widely used. 

During the last ten years, a number of unified constiutive models have 
been proposed which retain the inelastic strain as a unified quantity with- 
out aritifical separation into plasticity and creep components. These in- 
clude the models developed by Bodner (Refs. 11-13), Stouffer (Refs. 14-15), 

Krieg (Ref. 16), Miller (Ref. 17), Walker (Refs. 18-19), Valanis (Refs. 20- 
21), Krempl (Ref. 22), Cernocky (Ref. 23-24), Hart (Ref. 25), Chaboche (Ref.' 

26), Robinson (Ref. 27), Kocks (Ref. 28), and Cescotto and Leckie (Ref. 29). 

The applicability of these viscoplastic constitutive theories (mostly to 
high temperature applications) has been investigated by several researchers. 

Walker (Ref. 19) compared the predictive capability of several models (Walker, 

Miller and Krieg) for Hastelloy-X at 1800°F. More recently, Milly and Allen 
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(Ref. 30) provided a qualitative as well as quantitative comparison of the 
models developed by Bodner, Krieg, Walker and Krempl for IN100. Both Refs. 

19 and 30 conclude that these models generally provide adequate results 
for elevated isothermal conditions, they provide poor and overly-square 
results at low temperature, the material constants are often difficult to 
obtain experimentally, the resulting rate equations are "stiff" and sensi- 
tive to numerical integration, and the models do not provide any satisfac- 
tory transient temperature capability. Beek, Allen, .and Hilly (Ref. 31) 
have shown chat all the unified viscoplastic models mentioned above can be 
cast. in a functionally similar form (in terms of internal state variables). 

None of the published literature provides a thorough evaluation of cur- 
rent viscoplastic constitutive models with comparison to experimental re- 
sponse for complex input histories. Such an evaluation is difficult at pre- 
sent for many reasons, namely: 1) Material constants for most models are 

usually available only for a single material and often for a single temper- 
ature; 2) The experimental procedures given by model developers for deter- 
mining material constants from experimental data are often sketchy at best; 

3) Material constants for some models are often obtained by trial-and-error 
and are not based on experiments; and 4) There is a lack of good experimental 
data against which the models can be evaluated (thac is, test data which is 
significantly different from that used to generate the material constants). 

The purpose of the present paper is to report some preliminary evalu- 
ations of several of the unified viscoplastic models (Bodner, Krieg, Miller, 
and Walker). These four models are evaluated with regard to 1) their sen- j 

sitivity to numerical integration and 2) their sensitivity to slight changes 

in input material constants. j 

I 

i 

CONSTITUTIVE MODELS CONSIDERED 


The constitutive theories which have been studied to date include Bodner 's 
(Refs. 11-15), Krieg's (Ref. 16), Miller’s (Ref. 17), and Walker's (Refs. 
18-19). These particular models were selected for this initial study be- 
cause material constants for Hastelloy-X were available for three of the 
models. Other models are currently being considered as material constants 
become available. Each model is listed below in uniaxial form using a con- 
sistent notation as presented by Beek, Allen and Milly (Ref. 31). In Ref. 

31, it is shown that all of the current viscoplastio models considered may 
be written in uniaxial form as 

a =» E(e - a, - e T ) (1) 

where a is stress,. E is Young’s modulus, e is strain, is the inelastic 
strain (internal state variable), and is the thermal strain. Each vis- 
coplastic theory postulates a particular growth law for the internal state 
variable (s) and the inelastic strain is obtained by time integration of 
the growth law for a^, i.e. 
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a 1 = j a^Ct'ldt' 


( 2 ) 


where 


a 


dOj 

1 " dt 


a- L (e, T, a 2 . ot^. 


• • O 

m 


(3) 


In equations (2) and (3), c is time, T is temperature, ct 2 is the back stress 
(related to the dislocation arrangement and produces kinematic hardening 
or the Bauschinger effect), and is the drag stress (which represents 
the dislocation density and produces isotropic hardening). 

Bodner’s Theory 

The growth law for the inelastic strain in Bodner's model may be writ- 
ten in uniaxial form as 


a. - 


2 

ST L ° 


sgn(a) 


(4) 


where 


a, = m(Z, - a.)W - AZ 
0 1 3 p 


t (^) C 


(5) 


w = a a. 

P 1 


( 6 ) 


The quantities E, D 0 , n, m, Z p A, and r are material constants. As 
noted before, the variable CI 3 is similar to the drag stress used in 
many models (a measure of isotropic hardening or dislocation density) . 

It is noted that the model contains no parameter representing the back 
stress and cannot account for the Bauschinger effect in kinematic harden- 
ing materials. The material constants are tabulated for IN1Q0 at 1350®F 
(732°C) in Table 1 (taken from Ref. 1A). 
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Hiller's theory contains nine constants which are tabulated for Hastelloy- 
X at 1800°F (982°C) in Table 3 (see Ref. 19). 
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Walker's Theory 


Walker's nonlinear viscoplascic cheory can be cast in the following 
uniaxial form 


i° - a ?iV 


a. 


a- 


sg n(0 - a 2 ) 


(13) 



/n,R \" 

{i“ii k 

(n 3 + n^R) 2,n ^ 1+n ^ R ■+ ij 


+ n ? ja 2 - a 2 , 
o 


,m-l 


(14) 


a 3 = n 8 l a l 


n 9 i a il a 3 ~ n io ^ a 3 " a 3 ) 


(15) 


where R is the cumulative inelastic strain 

1 1 3a, 


R - J lit 1 " 


dt' 


(16) 


The general model requires sixteen constants (E, n, m, q, nj, 03, . . . , 
n , a 2 and a 3 (t=0). In determining the constants for Hastelloy-X at 1800°F 

(982°C), Walker made several simplifying assumptions [including 013 => cons- 
tant «* a3(t=0)] which reduces the number of parameters to those shown in Table 4 
(see Ref. 19). Further, the constants reported in Ref. 19 were developed from tests 
using strain rates in the range 10 3 to 10 6 sec 1 and strain ranges of iO.6%. 


NUMERICAL TIME INTEGRATION STUDY 


The integration of the constitutive relationship given by equations 
(1). (2) and (3) forms an integral and extremely important part in any nu- 
merical solution of a nonlinear field problem. It has been observed by 
many researchers that the coupled system of ordinary differential equations 
defining the state variables may be locally "stiff" and thus are sensitive 
to the time step size and numerical algorithm. The accurate integration 
of these stiff equations can be accomplished by various means: use of small 

time steps, higher-order or multi-point integration schemes, subincrementation 
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procedures (Refs. 33-35), "smart" algorithms which attempt to select appro- 
priate time steps in order to achieve accuracy and stability (Refs. 36,37), 
algorithms tailored for individual constitutive theories (Refs. 32,37), or 
combinations of these approaches. In general, the computation time required 
for the accurate solution of materially nonlinear problems is directly re- 
lated to the numerical integration scheme used. 

Regarding the constitutive models reviewed herein. Walker (Ref. 32) 
uses a stable, iterative implicit scheme which takes advantage of the func- 
tional ’form of the integrand in the development of the reccurence relation. 
Miller originally used Gear's method (Ref. 36) to integrate the stiff equa- 
tions in his theory but later concluded in Ref. 37 that an implicit back- 
ward difference method was more economical and preferable to either Gear's 
method or the explicit Euler forward integration method. The type of num- 
erical integration scheme used by Bodner and Krieg is not known. 

The selection of an appropriate time integration scheme to be used in 
a computer code is very important but is often based on the answers to such 
questions as: "What is available in the present code?", "What will work 

most of the time?", "What can we use that most users will understand?", 

"What is the cheapest and easiest to use?", and the like. The usual re- 
sponse given is "it depends on the problem being solved!" 

In general, equation (3) may be integrated between time t and t + At 
by writing 




dt 


(17) 


or 



% 


Acij =* cij(t + At) - a^(t) 


t+At 

I * 


t 


dt 


( 18 ) 


where is defined by the particular constitutive theory being used. The 
present investigation considers four integration schemes: explicit Euler 

forward integration, implicit trapezoidal method, trapezoidal predictor- 
corrector (iterative) method, and Runge-Kutta 4th order method. The approx- 
imations for each of these methods is given in Table 5. 

Each of the integration schemes in Table 5 were used to obtain stress- 
time and stress-strain responses for the four constitutive models considered 
herein when subjected to the uniaxial, alternating square-wave strain- rate 
history shown in Fig. 1. Figure I shows the 35 second response obtained 
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by Krieg's theory for Hastelloy-X at 1800°F using a time step of 0.1 sec- 
onds. For this time step, the Euler and trapezoidal predictor-corrector 
methods provide essentially the same results and are virtually identical 
to that obtained for all methods using a time step of 0.005 seconds. The 
4th order Runge-Kutta method generally overestimates the peak response while 
the trapezoidal method underestimates the response. Figure 2 presents re- 
sults for three integration methods such that the total computation time 
for a 35 second response solution is approximately the same. For equiva- 
lent computation times, the Euler method provides the most accurate results 
although smaller time steps are required. Similar results are observed for 
Miller's model. 

Figures 3 and 4 illustrate that various constitutive models may behave 
appreciably different using the same integration method (in this case the 
Euler method). In Fig. 3, Miller's theory (for Hastelloy-X at 1800°F) gives 
considerable oscillatory response for a time step of 0.005 seconds while 
Walker's theory shown in Fig. 4 gives a much smoother response for the same 
time step. Comparing Figs. 3 and 4, it is seen that a smaller time step 
is required (with Euler integration) in Miller's theory than in Walker's 
theory. 

Figure 5 presents results for IN100 at 1350°F using Bodner's model. 

Time steps were chosen for each integration scheme to obtain solutions which 
required approximately equal computation times. These results, when com- 
pared to solutions with much smaller time steps, indicate that the Euler 
method provides the most accurate results. Again, the time step used Is 
smaller than that for the other methods but the computation time is the 
same (for integrating the constitutive equations). 


SENSITIVITY STUDY FOR MATERIAL CONSTANTS 


In the previous section, results were presented which showed how the 
numerical integration method used to integrate the constitutive equations 
could affect the accuracy and computation times of predicted results for 
stress-time and stress-strain responses. In this section, we consider 
another important parameter in the application of any constitutive theory. 
Namely, "how does the accuracy to which material constants are determined 
from experimental test data affect the predicted response?" 

Figures 6 and 7 present results for Walker's model (Hastelloy-X at 
1800°F subjected to an alternating square-wave strain-rate history as shown) 
wherein specified input material constants have been adjusted by 5%. Fig- 
ure 6 shows the effect of a -5% change (error) in the stress exponent n 
(the most sensitive parameter). Figure 7 shows that a +5% error in all test 
data required to compute material constants results in significant predicted 
response errors, up to 30% over-prediction in the stress at a time of 35 
seconds (during the relaxation period) . 

Figures 8 and 9 present similar results for Krieg's model (Hastelloy-X 
at 1800°F) and Bodner's model (IN10Q at 1350°F), respectively. Both results 
indicate that the most sensitive parameter is the stress exponent "n" and 
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that a 5% error in specifying n may produce significant errors in the pre- 
dicted response. Miller's model appears to be much less sensitive to er- 
rors in input material parameters. 

Figure 10 provides a comparison of the Miller, Krieg, and Walker models 
for the Hastelloy-X test at 1800°F (using constants obtained by Walker for 
all models). The Euler method was used with a time step of 0.0005 seconds 
which provides a solution with no significant truncation error. The results 
obtained here show approximately 10-15% differences in peak stress ampli- 
tudes between the three constitutive models. Since no experimental results 
are available at this time, no conclusions can be drawn as to which model 
more accurately represents observed test data. However, the results do 
point out that significant differences (greater than 15%) can be obtained 
for stress peaks and stress relation values through the use of different 
constitutive models. 


CONCLUSIONS AND FUTURE WORK 


The results of this study are not complete since only a portion of the 
available constitutive models and numerical integration schemes have been 
considered. However, some tentative conclusions can be reached. First, 
it appears clear from the present investigation, and the work of others, 
that simple integration schemes (like the Euler forward defference method) 
are often preferable to more complex schemes from the standpoint of accur- 
acy, computation time, and ease of implementation. Although not reported 
herein, our work in progress indicates that Euler's method used with a simple 
subincrementation strategy provides the most accurate and economical solu- 
tion for most constitutive models. 

The sensitivity study on material constants indicates that most visco- 
plastic constitutive models are significantly sensitive to one or more mat- 
erial constants derived from laboratory tests. It has been shown that a 
5% "error" in laboratory measurements may lead to errors of 25%, or greater, 
in predicted stress responses. Although most model developers have fine- 
tuned their models and input material constants for specific material/temp- 
erature/strain-rate combinations, it is not clear that end-users will be 
able to do so when called upon to develop material constants for a new sit- 
uation. The problem can be negated to some extent by defining more explicit 
testing procedures for obtaining material constants and by guidelines de- 
fining which constants are most sensitive to experimental error. 

Our current and future work concerns the application of several inte- 
gration schemes to the other constitutive theories, investigation of sub- 
incremental strategies, and consideration of "smart" integration methods 
which detect local "stiffness" and adjust time steps but without signifi- 
cant computational expense. The material parameter sensitivity study will 
be continued by considering other constitutive theories, and more impor- 
tantly, by comparison with laboratory tests which involve complex thermo- 
mechanical loadings including transient temperature inputs. 
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Table 1. Material Constants Used in Bodner 's Model 
for IN100 at 1350°F (732°C) 


Bodners notation 

Beek and Allen's notation Numerical Value 

E 

E 

21.3x10* 

psi 

n 

n 

0.7 


Z 1 

m 

Z 1 

m 

1.105x10 
2. 57x10 3 

6 psi 
psi 1 

D o 

A 

D o 

A 

10 4 sec 
1 . 9xl0~ 3 

l 

sec 1 

r 

r 

2.66 



Z I 

0.6xl0 b 

psi 

e (t-0) 

a ]> (t=0) 

0.0 


Z o 

a 3 (t=0) 

0.915x10 

6 psi 


1 79 



Table 2. Material Constants Used in Krieg's Model 
for Hastelloy-X at 1800°F (982°C) 


Walker's notation Beek and Allen's notation Numerical Value 

for Krieg's constants 



Ci 

1.0 

n 

C 2 

4.49 

A-i 

a 2 

a 3 

C 3 

C 4 

Cs 

i.OxiO 6 psi 
6 . 21 x 10 6 psi 1 sec 1 
4. 027x1 0~ 7 psi -2 
100 psi sec 1/,n , 

4.365 psi 1_n sec l/n “ 2 

A 4 

As 

Ct 

c? 

E 

E 

13. 2xl0 6 psi 

Ko 

a3 o 

59,292 psi sec l/a 

n 

n 

4.49 

c(t» 0 ) 

ai(t» 0 ) 

0.0 

fl( t= 0 ) 

ct 2 (t= 0 ) 

0.0 

K(t-O) 

Cl3(t-0) 

59,292 psi 


Table 3. 

Material Constants Used in Miller 1 
for Hastelloy-X at 1800°F (982°C) 

1 s Model 

Miller's notation 

Beek and Allen's notation 

Numerical Value 


n 

n 

2.363 

B0' 

B0' 

2.616x10 6 sec 1 

Hi 

Hi 

1 x 10 s psi 

Ai 

Ai 

1.4053xl0~ 3 psi " 1 

h 2 

h 2 , 

100 psi sec 1 ' 11 

c 2 

c 2 

5,000 psi 

a 2 

a 2 

4.355x10 12 psi 3 

E 

E 

13.2x10 s psi 

e(t-O) 

01 (t= 0 ) 

0.0 

R (t= 0 ) 

a 2 (t= 0 ) 

0.0 

Do 

0.3 (t* 0 ) 

8,642 psi 
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Table 4. Material Constants Used in Walker's Model 
for Hastelloy-X at L800°F (982°C) 


Walker's notation 

Beek and Allen's notation 

Numerical Value 

o 

Q 

ct 2 

-1,200 psi 

ni 

ni 

0 psi (not used) 

n 2 

n 2 

1x10 s psi 

n 9 

* 

312.5 

n? 

«7 

2.73x10 3 psi 1 m sec 

n 

n 

4.49 

m 

at 

1.16 

E 

E 

13.2x10 s psi 

c(t-O) 

ai(t=0) 

0.0 

Q(t=0) 

aj(t-O) 

0.0 

K(t=0) 

a 3 (t»0) 

59,292 psi 


ne,ng,nio,q 

0 (not used) 


* 


3 

SR. 


(n 3 + n^R) 


In 


/ n 5 a 

\i« 6 r 



Table 5. Numerical Integration Approximation for Aci.^ « a 1 dt 


Method 


Approximation 

Euler Forward Difference 

Ac^ ° At (t) 

Trapezoidal 

Rule 

AOj," (t) + (t+At) ] 

Trapezoidal 

Predictor-Corrector 

Same as trapezoidal except iterate 

Runge-Kucta 

4th Order 

Aa x - fa + 2K 2 + 2K 3 + K 4 ) 



1^ = At a 1 (t,a 1 (t)) 



K 2 =» At a 1 (t+At/2,a 1 (t)+K 1 /2) 



= At a 1 (t+At/2,a 1 (t)+K 2 /2) 



= At (t+At.a^ (t)+K 3 ) 
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Fig. 1 Comparison of integration methods for Krieg's theory (Hastelloy-X 
at 1800°F) 
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•Stress (KSI) 



g. Stability and accuracy of Fig. 5 Comparison of integration methods 

Euler s method for Walker's for Bodner's theory with equal 

theory (Hastelloy-X at 1800 F) computation time allowed for each 

method (IN100 at 1350°F) 
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Fig. 6 Sensitivity of Walker's theory to -5% change in input constants 
(Hastelloy-X at 1800°F) 
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Fig. 10 Comparison of stress-time predictions for Miller, Krieg and 
Walker theories (Hastelloy-X at 1800°F) 
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ABSTRACT 


The sensitivity of the unified constitutive theories, proposed by 
Bodner, Walker, Krieg, and Miller, to numerical integration techniques 
and slight changes in material parameters was investigated. Evaluations 
were based upon numerical simulations of Hastelloy-X at 1800° F in uniaxial 
form by specifying input strain histories and comparing output stress 
histories. The integration sensitivity study involved the use of both 
single and multi-step integration schemes. The various algorithms were 
compared with regard to accuracy, stability, and computational economy. 

The material parameter sensitivity was studied by varying the material 
constants by a specified amount and comparing predicted responses. 

Numerical comparisons show that, of the numerical integration methods 
studied, Eulers method is the most accurate, stable, and efficient pro- 
cedure. The input sensitivity studies indicated that some constitutive 
models are more sensitive to experimental errors than others and that a 
5 l error in certain material constrants may lead to 15-30% changes in 
predicted response. 
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CHAPTER I 
INTRODUCTION 

For many years there has been a substantial amount of research 
performed in the area of modeling the constitutive behavior of in- 
elastic deformation. Recently, new unified formulations have been 
proposed which are attracting much attention because these theories 
treat the inelastic strain as a unified quantity, incapable of being 
separated into time dependent and time independent parts. However, 
these new models have two potential problems which can effect their 
predictive capabilities. These include sensitivity to the numerical 
integration of the differential equations describing the inelastic 
strain component, as well as a difficult and sometimes ill-defined 
method of determining the material parameters. 

Recently, Beek, Allen and Milly [1] have shown that a number of 
proposed thermoviscoplastic constitutive models can be cast into a 
common internal state variable framework. The first state variable 
is the inelastic strain (e^) whose rate of growth is a function of both 
current conditions (stress o and temperature T) and one or more internal 
state variables. These unified formulations form a system of nonlinear 
ordinary differential equations in the independent variable time (t). 
Many researchers [2-16] have observed these equations to be well 
behaved in some portions of the strain response and mathematically 
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stiff in other regions. Therefore, in order to obtain a stable and 
accurate prediction of inelastic behavior, one must carefully consider 
the method of integration. There are a variety of numerical 
algorithms available, from the simple Euler integration scheme 
wherein one must judiciously select an appropriate time step size to 
insure solution stability and accuracy, to the more sophisticated 
method of Gear [7-9] where the step size is automatically selected 
after checking for the existence of stiffness. 

The growth laws which are used to describe the stress-strain 
response of a particular material rely on a number of material 
parameters. The material parameters (or constants) are determined, 
either implicitly or explicitly, from a series of uniaxial tests 
which are prescribed by the model developer. Typically some 
combination of the following tests is required: 1) constant strain 

rate tensile tests; 2) constant strain rate tensile tests with 
intermittent hold times; 3) constant strain rate cyclic tests; 

4) relaxations tests; 5) creep tests; or 6) stress drop tests. 

The application of test data in determining material constants varies 
from model to model. In most cases, because of the lack of 
standardization, it would appear that the end user will be burdened 
with a trial -and-error method for determining specific material 
parameters. This burden might be eased if. a quantitative assessment 
of the model sensitivity to slight changes in material constants was 
provided. However, there is no such information available in the 
literature, nor is the problem even addressed. 
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The purpose of this research will be to evaluate numerical 
simulations of several unified constitutive models, in uniaxial form, 
by specifying input strain histories and comparing output stress 
histories. Integration methods and slight changes to material 
parameters will be the prime variables. This will allow for a 
qualitative analysis of solution stability, accuracy, and sensitivity. 

In fulfillment of these objectives, this thesis will proceed 
in the following manner. First, a review of the pertinent literature 
is presented in addition to a brief overview of the origin of 
stiffness. Next, the constitutive theories considered herein are 
discussed and the various methods used to determine their material 
parameters are outlined. Finally, the numerical integration and 
parameter sensitivity study is presented, wherein the use of 
various integration algorithms selected for this study are 
investigated and the effects of varying material constants, 
simulating "experimental error", are analyzed. Evaluations are 
based upon numerical simulations of Hastelloy-X at 1800° F, in 
uniaxial form, when subjected to several strain input histories. 
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CHAPTER II 
BACKGROUND 


Constitutive Modeling 

Constitutive modeling of inelastic deformation is not a new 
concept, it has been in practice since 1864 starting with the 
early work of Tresca [17]. Levy [18] and von Mises [19] added to 
the study of plasticity by developing what have come to be known as 
the Levy-Mises equations. The development of classical plasticity 
continued with the work of Hencky [20], Prandtl [21], Prager [22], 
and Ziegler [23] to name just a few. These early models were 
incremental in nature, assumed that plasticity and creep could be 
separated, and they incorporated a yield function, flow rule, and 
hardening rule to define the plastic strain increment. 

In addition to classical plasticity based-theories, other 
fundamental bases were proposed to describe inelastic deformation. 
These include the microphenomenologically-based theories and thermo- 
dynamical ly-based nonlinear viscoelastic theories. The aforementioned 
models can be further categorized as unified or uncoupled theories, 
where the two differ in their treatment of rate-independent and 
rate-dependent strain components. The uncoupled theories partition 
the inelastic strain into plastic and creep components; whereas, 
the unified constitutive models retain the inelastic strain as a 
unified quantity. However, it should be noted that the partitioning 
method used in the uncoupled models is questioned because it has no 
formal micromechanical basis and does not account for creep and 
plasticity interaction. 
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Over the past decade a number of unified constitutive theories 
have been developed. These new theories include the work of Bodner, 
et al. [24-26], Cernocky and Krempl [27,28], Cescotto [29], Chaboche 
and Cailletaud [30], Hart [31], Krempl [32], Krieg, et al. [33], 
Kocks [34], Miller [35], Robinson [36], Stouffer [37,38], 

Valanis [37,40], and Walker [41,42]. 

As stated earlier Beek, Allen, and Milly [1] and Krieg [16] 
have shown that all thermoviscoplastic constitutive models can be 
cast into a common internal state variable framework. For the 
uniaxial case, a general form of these unified models may be 
written as: 


a 

= E(e - e* - e^) 

(1) 

•I 

e 

- f ( e , T , B , D , , . . , a^) 

(2) 

B 

= hg e 1 - r B (T,B) 

(3) 

• 

D 

= hg e 1 - r 0 (T,D) 

(4) 


where o is the applied stress, E is Young's modulus, e is the strain, 
e * is the inelastic strain or first state variable (with both 
plasticity and creep included), eJ is the thermal strain, B is the 
back or rest stress corresponding to the dislocation arrangement 
which produces kinematic hardening or the Bauschinger effect, and 
D is the drag stress corresponding to the dislocation density which 
produces isotropic hardening. In addition, f is the inelastic strain 
rate function, T is the temperature, c« are additional state variables 
which are presently undefined, and hg, hp, rg, and r^ are the hardening 
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and recovery functions for the back or drag stress (designated by the 

subscript B or D), respectively. A superposed dot above the 

variables denotes differentiation with respect to time. 

The abstract model presented above (equations (l)-(4)) contains 

all three of the commonly used internal state variables; however, 

specific theories may omit one or more of the growth laws which 

characterize either the back stress, drag stress, or other variables. 

In all cases, the inelastic strain is obtained via time integration 

of the specified internal state variable growth law for e*, i.e., 

t 

^(t'Jdf , (5) 

where e 1 is defined in functional form by (2), t is the current time 
of interest, and t' is a dummy variable of integration. 

Integration Methods 

A number of integration strategies have been evaluated 
[2,3,5,6,10-14] for integrating constitutive equations which have 
"locally" stiff regimes. These algorithms have been implemented 
into both large multi-axial finite element programs and small 
uniaxial constitutive codes. For the purpose of this discussion, 
only those publ ications dealing with integration on the constitutive 
level will be reviewed. 

Numerical schemes for the solution of the inelastic response 
behavior were explored by William [2]. He used a generic constitutive 
formulation in order to examine the algorithmic properties of various 
numerical solution methods. These included the explicit forward 
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Euler method, the forward gradient method, the predictor-corrector 
method, and the Newton-Raphson method. A qualitative comparison of 
the four incremental solution methods was presented for pure creep 
and pure relaxation behavior; however, no specific conclusions were 
stated with regard to the "best" scheme. 

An explicit trapezoidal method for handling deformation in 
stiff regimes was proposed by Shih, Delorenzi, and Miller [3], 

Results using the explicit trapezoidal scheme demonstrated both 
improved stability and computational efficiency over either explicit 
forward Euler integration or the trapezoidal method with a Newton- 
Raphson corrector. It was shown that the predictor-corrector 
method eventually converged to the same solution predicted by the 
explicit trapezoidal scheme, while the Euler predicted response was 
unstable. 

Another promising numerical technique, called the a or NONSS 
(Noniterative, Self-correcting Solution) method, has been developed 
by Tanaka and Miller and is discussed in Reference [4]. The implicit 
quantities are removed by a Taylor series expansion of the prescribed 
growth laws through the use of an integration operator a. The 
integration operator, a, has the range of (0.1). This alogrithm 
also incorporates automatic time-step control and an error corrector. 
Another scheme similar in style to the aforementioned is the ©-method 
(described in Reference [5]) which incorporates the explicit forward 
Euler method (©=0), the implicit trapezoidal scheme (0=^), or the 
implicit Euler strategy (©=1). Numerical experiments indicated that 
the forward Euler integration scheme (©=0) with an automatic stepping 
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and error control was computationally more efficient. 

For integrating constitutive equations of the work hardening- 
recovery type, Miller and Shi h [6] have developed a method wherein 
improved accuracy or larger time steps or some combination of both 
can be obtained. This is accomplished through a Taylor series 
expansion of the recovery function resulting in a more accurate 
calculation of the internal variable(s) during each time step. 

A comparison between this special algorithm and Euler integration 
from Zircaloy-2 simulations verify the accuracy and stability of 
this method. 

Gear's method [7-9] has also been an effective tool for 
integrating stiff differential equations. It is a multi-step 
predictor corrector scheme whose order is automatically chosen. 
Either an Adams' method or methods suitable for stiff equations 
can be selected. Three alogrithms are available, STIFF-0, STIFF-1, 
and STIFF-2. Chang and Chang [5] believe that Gear's package is 
well suited to solving one-dimensional constitutive relations, but 
is much too cumbersome for use on the structural level. Based on 
this work, it may be tentatively concluded that Gear's method is not 
suitable for implementation into finite element packages. 

A very comprehensive study on numerical integration of stiff 
constitutive models was reported by Kumar, Morjaria, and 
Mukherjee [10]. Hart's equations for predicting inelastic behavior 
were used to compare various numerical approximation strategies. 
Several input histories of stress and strain were considered and 
the accuracy and computational efficiency of the results were 
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compared. The study concluded that a simple Euler type one-step 
method with automatic time-step control worked well; however, it 
emphasizes that a one-step strategy is not necessarily the best. 

A combination of the two-step Adams’ method (for outside the 
viscoplastic limit) and Euler's method (for inside the viscoplastic 
limit), was found to be the most efficient in terms of computational 
speed and accuracy. 

The Origin of Stiffness 

The purpose of this section is to familiarize the reader with 
the concept of stiff equations through a physical and mathematical 
definition. By understanding this concept, researchers may then 
devise methods to test for stiffness or more simply concede to its 
existence and develop appropriate solution algorithms. 

A system of differential equations is said to be stiff if the 
physical processes being modeled contain time constants with 
different scales. In the case of a uniaxial test specimen, there 
is grain boundary or dislocation movement in one time scale; 
whereas the total stress-strain response is measured using a 
different time constant. The stiffness of equations (2) thru (4) 
arises from the following phenomenaon. The inelastic strain rate, 
e*, may be a strong function of o, B, and D in (2); that is small 
changes in o, B, and D cause large changes in c*, which then 
influences a, B, and D through (1), (3), and (4). Ultimately, the 
absolute stability of equations (2) thru (4) requires that the step 
size used during numerical integration be not much greater than the 
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smallest time constant present in the physical system; however, for 

computational reasons this is not always possible. Therefore, 

integration methods must be used that are "stiffly stable". 

The mathematical definition as described in [44] is as follows: 

"The nonlinear system Y' = F(X, Y) is said to be stiff in an 
interval I if for X e I the eigenvalues A. of the Jacobian J 
satisfy the following conditions." K 


1 ) Re A^ < 0, k = 1 .... ,m 

2) max . | _ min )e . , . 

l<k<iJ iJ l<k<nJ 

Where the a, , k=l,...,m are the eigenvalues of F. 
The Ratio 


max I 
l<k<m' 

mTn 

l<k<nr 


Re A k j 
Re x k | 


is called the stiffness ratio 


( 6 ) 


Basically this means that if the real eigenvalues are negative or 
if the spread of real eigenvalues is large, then the equations are 
stiff. Gear's method [7-9] uses this definition to test for the 
existence of stiffness. However, this can be very expensive 
expecially when large systems of equations are being analyzed. 


Material Testing 

Each constitutive model has a number of material parameters which 
are determined through a set of complex tests; therefore, the end 
user must be informed about the effect that small variations, i.e. 
typical experimental error, have on the particular models' predictive 
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capabilities. In general, the testing of materials is a difficult 
process in which the propensity for error is great. This experimental 
error is known to occur in all stages of testing and will be discussed 
herein; however, the reasons for its occurrence are beyond the scope 
of this research. 

Experimental error may develop in all phases of testing; however, 
it occurs most frequently in the following areas. The nonrepeatability 
of strain measurements for uniaxial strain rate, cyclic or creep 
tests makes interpretation of results difficult even when new high 
technology axial and diametral extensometers are used. When testing 
requires elevated temperatures, error may result while measuring and 
maintaining the specimen's temperature. Finally, experimental error 
occurs in its most rudimentary form in the statistical variation of 
the test sample's measurements. While these are only an example of 
where experimental error occurs, they are representative of the 
problems experimentalists must cope with and why absolute values of 
material constants are difficult if not impossible to obtain. 

Objectives 

The objective of this research is to study the sensitivity of 
the unified constitutive theories proposed by 1} Bodner, et al., 

2) Krieg, Swearengen and Rohde, 3) Miller, and 4) Walker to 
numerical integration techniques and slight changes in material 
parameters. Evaluations are based upon numerical simulations of 
Hastelloy-X at 1800° F in uniaxial form by specifying input strain 
histories and comparing output stress histories. The integration 
sensitivity study involves the use of both single and multi step 
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integration schemes. The material parameter sensitivity is studied 
by varying material constants by a specified amount and comparing 
predicted responses. 
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CHAPTER III 

CONSTITUTIVE THEORIES CONSIDERED 

The constitutive models selected for this study include the work 
of Bodner (et al.) [24-26], Krieg (Swearengen and Rohde) [33], 

Miller [35], and Walker [41-42]. These models were chosen primarily 
because of the availability of material parameters for the same 
material (Hastelloy-X) , the exception being Bodner 's model (constants 
for Bodner's model were obtained by simulating the required uniaxial 
tests using Walker's model). A secondary reason for selecting these 
models is that they appear to be the most qualitatively attractive 

theories available [42, 43] and thus warrant further evaluation. 

In the following sections, each model is presented in one-dimensional 
form and the methods used to compute the material constants are 
discussed. 


Bodner, et al . 1 s Theory 

The model proposed by Bodner, et al . is a microphenomeno- 
logically based theory for characterizing inelastic behavior. This 
theory was developed to predict the response of a material which work 
hardens isotropically. The model contains two internal state variables 
the inelastic strain and the drag stress. Bodner's model contains 
no parameter which represents the back stress and thus cannot 
account for the Bauschinger effect in kinematic hardening materials. 

The uniaxial-differential form of Bodner's model may be written as 


•I 

e 


_2 

/3 


(-< 


n±l)(£) 2n> \ 

2n m D ; 


sgn ( 0 ) 


(7) 
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where 


D = 


m(D 1 


- D)W p - AD ] 


D-D, 




•I 

a c 


( 8 ) 


(9) 


The material parameters to be determined are E, e Q , n, m. A, r, 

D q , D^ , and D^, where D q is the initial value of D„ The constant e Q 

represents a limiting value of the inelastic strain rate and is usually 
4 -1 

assumed to be 10 sec unless the strain rates are very high. 


Determination of Material Parameters ., The material parameters 
for Bodner's model are determined through a series of constant strain 
rate tensile tests and constant stress creep tests,, Control 
parameters for both sets of experiments should encompass values 
which are compatible with the numerical simulation to be performed. 

The first step in the evaluation of the material constants, as 
described by Stouffer [37] and Milly and Allen [43], is to determine 
n from the constant strain rate tensile data. By observation, when 
the stress reaches a saturated value during a constant strain rate 
test equation (7) is satisfied if and only if D is a constant* 
Therefore, by neglecting recovery via rapid loading of the test 
specimen, D can be assumed to be in its fully work hardened state 
and have a maximum value D-j. Thus (7) may be rewritten, to account 
for the steady state flow condition, as 


In [-ln(^| 


■ I 


-)] = -(2n)lno 


+ [2n lnD^ln^)] 


( 10 ) 


If the experimental data base is good, the left hand side of equation 
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(10) must be a linear function of lna„ Thus, a plot of these two 
variables (as depicted in Fig* 1(a)) yields a straight line where 
the slope and intercept are defined by the quantities -2n and 
[2n In D.j + ln((n+l)/2n)], respectively, from which n andD^, may be 
determined. 

To determine the parameters m and D Q , constant strain rate 
tensile test data is used. If work hardening is assumed negligible, 
(8) may be rewritten as 


dD = m(D 1 -D)dW p , (11) 

which can be integrated to give 

ln(D r D) = -mW p + ln^-D^ , (12) 

where W p is the inelastic work defined by (9) and D q is the initial 
value of hardness D. To obtain D for measured values of stress and 
inelastic strain rate, equation (7) is inverted and may be written as 


D = a 


/ 2n , . /Ox 


l/2n 


(13) 


Therefore, by using values of o and e* from the nonlinear portion of 
the tensile response data the quantity D may be determined by 
equation (13). Since equation (12) is a linear representation, of 
ln(D-|-D) and W p , a plot of these variables will be a straight line 
(see Fig. 1(b)) having a slope of m. Thus, knowing m, D Q may be 
determined through simple substitution and rearrangement of equation 


( 12 ). 




In (o) 


a) Evaluation of Strain Rate b) Evaluation of Hardening 

Exponent n. Parameter rru 



c) Evaluation of Hardening 
Recovery Parameter r. 


Fig. 1. Plots Used to Obtain Material Parameters 
For Bodner ' s Model . 
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The last set of constants. A, r, and are computed using data 
obtained from constant stress creep tests. During secondary creep 
the creep rate is approximately constant; therefore, the value of D 
must be constant and D = 0. Thus, the hardening rate equation ( 8 ) 
for secondary creep becomes 

/D-D.\ r 

m(D r D)W p = AD, (- 5 ^) . (14) 

linearization of (14) results in 

In [m(D r D)W p ] = r ln(^) + lnfAD^ , (15) 

where stable values of D are computed using equation (13) and D£ is 

• 

the minimum calculated value of D. By plotting ln[m(D.-D)W ] as 
D-D p 

a function of ln( — ^ - ) , similar to Fig. 1(c), the variables A and r 

U 1 

may be determined. From equation (15), it may be seen that r 
represents the slope of the straight line plot and 1 n (AD-j ) i s the 
intercept. Hence, A may be computed once the value of the 
ln[m(D^-D)W p ] intercept is known. 

Since no experimental data for Hastelloy-X at 1800° F was 
available, the material parameters for Bodner's model were computed 
using Walker's theory (as described in the next section) to predict 
the deformation response. Numerical simulations included both 
constant strain rate and constant stress (creep) tests. The material 
constants computed for Bodner's model are tabulated in Table 1. 
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Table 1. Material Constants Used In 

Bodner's Model For Hastelloy-X 
at 1800° F 


Bodner's Constants 

Numerical Value 

E 

13.2 x 10 6 psi 

e 0 

1 .0 x 10^ sec" 1 

n 

0.103 

m 

- 0.178 psi 

r 

0.294 

A 

0.162 x 10" 1 sec" 1 

D 1 

0.329 x 10 7 psi 

D 2 

0.305 x 10 7 psi 

e’(t=0) 

0.0 

D(t--O) = D q 

0. 31 7- x 10 7 psi 
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Walker's iheory 


The model proposed by Walker is a viscoelastic theory developed 
by modifying the constitutive relation for a three parameter visco- 
elastic solid. This model contains three internal state variables; 
the inelastic strain, the back stress, and the drag stress. The 
back stress introduces nonlinear kinematic hardening into the model 
to account for the Bauschinger effect and the drag stress introduces 
isotropic hardening into the model to account for cyclic hardening 
or softening of the material. Creep, relaxation , and strain rate 
effects are modeled by a power law function for the inelastic strain 
rate. The uniaxial-differential form of Walker's model may be 
written as 

•I = 1 ° - B -L n 1 ( 0 - B) (16) 

D n 

B = (n ] + n 2 ) e 1 - (B - B - n ] t 1 ) G (17) 


D = n 8 | e 1 | - n g \z l \ D - n ]0 (D - D Q ) q , 


(18) 


where the value of G is defined to be 

A _ i • 1 1 3 


3R 


/ M 

(n 3 + n 4 R) In^^Tp + 1 


+ n 7 | B-B | m_ 1 (19) 


and R is the cumulative inelastic stran 



Walker made several simplifying assumptions to equations (17) and (18), 
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including D = constant for the Hastelloy-x material. Therefore, 
history dependence in this model is exhibited entirely by a 
simplified growth law for the back stress. 

B = (n 1 + n 2 ) e 1 - (B-B-n^e 1 ) (nge; 1 + n ? | B-B j m_1 ) (21) 


Thus, the material parameters to be determined are E, n, m, n-j , 

O 

n 2 , n 7> rig, B, and D. 


Determination of Material Parameters . The material parameters 
for Walker's model are determined through a series of constant strain 
rate steady state hysteresis loops under fully reversed strain 
controlled conditions at a constant temperature. 

The first step in the evaluation of the material constants, as 

O 

described by Walker [41], is to determine the value of B. The 
asymptotic steady state value of B may be obtained under conditions of 
uniaxial tension loading at a known constant strain rate because 
the back stress will saturate to a constant value as e 1 ” if n^ = 0. 


I 


li m 


n^ e 


g _ g _|_ 

£ I_ > - n Q R + n,lB-B| m ' 1 


( 22 ) 


The back stress becomes rate-independent and will saturate to a 
maximum value for rapid monotonic loading in tension and compression 

• . T 

when R = e . Thus, the static thermal recovery term involving n 7 
in (22) may be neglected and maximum values of the back stress, in 
tension and compression, may be obtained by 


max = B + V n 9 


(23) 
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and 


B = B - n 0 /n a 
max 2 9 


( 24 ) 


For large strains, stress is assumed to approach an asymptotic 
such that e* = e; therefore, equation (16) may be rewritten as 


a = B + D 


1/n 


(25) 


The maximum values of stress in rapid tension and compression loading 
are then obtained as 


“max ' B * n 2 /n 9 + 0 


(26) 


and 


“max ■ B - V"9 - D “ ,/n 


(27) 


Hence, B may be obtained using the relation 


n _ 1 / t c , 

B - 7T (a + a ) 
2 ' max max' 


(28) 


To determine the constants n and D equation (25) may be 

rewritten in terms of steady state values of o and B in tenson as 

. n 

e = \ k / . (29) 


t p t 
a - b 
max max 


D 


If o-j and 0 ^ denote the corresponding maximum tensi le stresses at the 
two maximum strain rates, ^ and where B attains its maximum 
saturated value, then equation (29) gives 

n - in (e,/; 2 )/1n j <V B L>' ( vC’ | (30) 



22 


and 

° * h - O v ,/n • < 3i > 

The value of B* ax used in equations ( 29 ) - ( 31 ) is determined 

from stress drop tests. If the material is cycled under steady 

state strain controlled conditions around a closed hysteresis 

loop as depicted in Fig. 2, positive creep will occur at point C 

if the stress is held constant and alternately if the stress 

is held constant at point I, negative creep will occur. Since the 

saturation stress in tension, B t , is desired, a point F on the 

max 

hysteresis loop between points C and I must be found such that the 
creep rate ceases instantaneously when the stress is held at a 
constant value. Walker indicates that it is difficult to estimate 
the value of B^ ax using the above procedure because the creep 
rates are usually small on the unloading branch. He suggests that 
a better estimate may be obtained by a relaxation test, noting 
the point at which initial relaxation changes from positive (stress 
decreases) to negative (stress increases), since relaxation 
progresses more rapidly than creep. 

Values for the material parameters n^ and ng may be determined 
from the initial monotonic stress-strain curve. For rapid loading 
rates n-| = 0 and G = ng R, therefore the back stress may be expressed 
as 

B = B + (B 1 - B) (1 - e“ ' o/E ^) 

v max ' v ' 


( 32 ) 
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For tensile loading equation (16) may be written as 


£ 

E 


1_ d£ • 
E de 


f -i 1 dOx _ /0“Dx 

£“i- £“rj £ “ Y (1 - £ jjg.) ( n ) 


o^JB\ n 
D 


(33) 


and combined with equation (32) to give 


n g = (e - o/E) 


-1 


In < 


B L - B 
max 


<4x ' « + D * 1/n <’ 


1 da xl/n 

E de ' 


r (34) 


where y is the applied strain rate and da/de denotes the instantaneous 
slope of the uniaxial monotonic stress-strain curve at the point 
where the stress is a and the strain is e. An estimate of the 
constant can then be made by rewriting equation (23) as 


n 

2 


(B 


max 


- B) 


(35) 


The final two constants to be determined are m and n^. 
Equation (21), with n^ = 0, yields the expression for the back 
stress rate for tensile loading as 


B = n 2 e 1 - (B - B)(n g e 1 + n^B-B)" 1 ' 1 ) . 


(36) 


At the peaks of the hysteresis loops B saturates to a constant 
value, _B = 0 and o = 0 for tensile loading, so that equation (36) 
gives 

n 2 e 4 = (B 4 - B)n g e 4 + n 7 (B 4 - B>" (37) 


or 


n 2 ^5 = ^ B 5 ‘ B ) n g ^5 + n 7^ B 5 “ B ) 


m 


(38) 
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If a. and o,- denote the maximum tensile stresses at the two 
smallest strain rates and then and B^, which represent 
the saturated back stresses for these two strain rates, may be 
determined by equation (25). By combining equations (37) and (38) 
m is obtained in the form 


m 


In 


n 2 £ 4 
. n 2 e 5 


In 


“ n g e 4 


n 9 e 5 






(39) 


Knowing m, equation (37) may be solved for n^. Note that n^ is 
assumed to be zero in these applications. 

The material constants computed for Walker’s model, simulating 
Hastelloy-X, were evaluated by Walker [41] and are tabulated in 
Table 2. 


Krieg, Swearengen, and Rohde's Theory 

The model proposed by Krieg, Swearengen, and Rohde uses a 
power law function to model the inelastic strain rate and is based 
upon the microstructural aspects of deformation mechanics. This 
model contains three internal state variables; the inelastic strain, 
the back stress, and the drag stress. The growth laws for the 
back stress and drag stress are formulated as the differences in 
hardening and recovery rates. The back stress is a kinematic 
hardening state variable which is used to model the Bauschinger 
effect and the drag stress is an isotropic hardening state variable 
which attempts to account for the dislocation density or mechanical 
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Table 2. Material Constants Used In Walker's 
Model For Hastelloy-X at 1800°F 


Walker's Constants 

Numerical Value 

E 

13.2 x 10 6 psi 

n 

4.49 

m 

1.16 

n l 

0.0 psi 

n 2 

1.0 x 10 6 psi 

n 7 

2.73 x 10' 3 psi 1 '" 1 sec" 1 

n 9 

312.5 

O 


B 

- 1 ,200 psi 

D 

59,292 psi 

e 1 (t=0) 

0.0 ■ 

B(t=0) 

0.0 

D(t=0) 

59,292 psi 
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strength of the material. The uniaxial-differential form of 
Krieg's model may be written as 


• I = 1< L_ - B] P (o . B) { 40) 

D n 

O 

° ° - 2 

8 = A, c - A 2 (B - 8) |B - B|(e A 3 (B_B) -1) (41) 

D = A 4 |e'| - A 5 (D - D 0 ) n . (42) 


It should be noted that equations (16) thru (18) vary in form from 
those originally proposed by Krieg [33], instead the notation 
adopted by Walker [42] is presented. The material parameters to 

o 

be determined are E, n, A 1 , A^, A^, A 5 , B, and D Q . 

Determination of Material Parameters . The material parameters 
for Krieg's model are determined through a series of constant 
strain rate cyclic tests and stress drop tests. 

o 

The constants n, Dq, B, and A-j in Krieg's model are analogous 

o 

to the material parameters n, D, B, and (n-j + n^) respectively, 
in Walker's Theory. Therefore, the interested reader is referred 
back to Walker's model for a detailed explanation of their 
determi nation . 

Under steady state conditions where o = 0, equation (40) shows 
that B = constant for constant strain rates. That is, 

8 1 = a] - D c] /n (43) 

or 

r - n ‘ Vn 
“2 - °2 ~ ^ e 2 


( 44 ) 
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Where arid are saturated values of the back stress for strain 
rates and respecti vely. Therefore, by equating equation (41) 
to zero and noting that e = e* if a = 0, gives (for c positive) 



= (e A 3 (B - B) 2 .,, 

t 


(45) 


By applying known values of the back stress (B^ and B 2 ) for the 
given strain rates (e^ and 1^) to equation (45) to give two 
equations with the unknowns, (A-j/A^) and A^, the constant (A^/A^) 
may be eliminated so that A^ is determined by 


Ini 


1 + 


A 3 


B 2 - B 
'2 \ B i " B 


(eW - B > 2 -l) 


(46) 


(B 1 - BY 


which can be solved using an interative technique. 

The constant A^ may be determined from steady state hysteresis 
loops where B = 0. Under these conditions, equation (41) may be 
written as 


A 2 " 


1 e 


(B 


- B) 2 (e A 3 (B - B > -1) . 


(47) 


The final two constants to be determined are A^ and A,.. The 
parameter A^ controls the amount of isotropic hardening produced 
by a given amount of strain. As such, it plays an important role 
in three different transient situations: 1) When a tensile 



29 


test is simulated, controls the rate of work hardening; 

2) When a creep test is simulated, controls the amount of 
primary creep; and 3) When a fatigue test is simulated, A^ controls 
the rate of cyclic hardening (or softening). A best-fit value for 
A^ can be tried and the resulting simulations compared with the 
test data for the same loading, to select the one which causes 
the closest agreement. Once a value for A^ has been selected, the 
constant A^ may be determined from the following relation 


fl 5 = 


V 


< D r D o>" 


(48) 


where is the steady state value of the drag stress at a constant 
strain rate, e c 

The material constants computed for Krieg's model, simulating 
Hastelloy-X at 1800° F, were evaluated by Walker [42] and are 
tabulated in Table 3. 

Miller's Theory 

The model proposed by Miller [35] is a phenomenologically based 
theory which uses a hyperbolic sine function for the inelastic 
strain rate. This model contains three internal state variables; 
the inelastic strain, the back stress, and the drag stress. The 
theory is capable of simulating many of the structurally significant 
deformation phenomena including cyclic hardening/softening, 
Bauschinger effect, strain rate effects, and annealing, among others 
The growth laws for the back stress and drag stress are formulated 
as the difference between hardening and recovery rates. 
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Table 3. 

Material Constants Used In Krieg's 
Model For Hastelloy-X at 1800° F 

Krieg's Constants 

Numerical Value 

E 

13.2 x 10 6 psi 

n 

4.49 

A 1 

1.0 x 10 6 psi 

CJ 

c 

6.21 x 10” 6 psi” 1 sec” 1 

A 3 

4.027 x 10” 7 psi" 2 

fl 4 

100 psi sec 1/,n 

A 5 

4.365 psi 1 " 0 sec 1/n ” 2 

B 

- 1 ,200 psi 

D o 

59,292 psi sec 1/n 

e 1 (t=0) 

0.0 

B(t=0) 

0.0 

D(t=0) 

59,292 psi 
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The uniaxial-differential form of Miller's model may be written as 


e = B^e' < sinh 


a-B 


VI. 5 1 


sgn(o-B) 


(49) 


(50) 


§ = H-j e 1 - H 1 B c o t { sinh A ] |B| } n sgn (B) 

D = H 2 1 e 1 1 (C 2 + j B j - A 2 D 3 /A 1 ) - H 2 C 2 B c e' {sinh(A 2 D 3 ) } n . (51) 


The material parameters to be determined are E, B^q' , n, , H 2 , 
A 1 ’ A 2’ S’ anc * D 0* 0n ^ cons ^ ants which depend on temperature 

are the coefficients B^G' and Dq. 


Determination of Material Parameters. The material parameters 
for Miller's model are determined through a series of constant strain 
rate steady state hysteresis loops under fully reversed strain 
controlled conditions at a constant temperature. 

The starting point for the calculation of the constants is to 
find the values of B^g', n, and A that satisfy one following steady 
state equations 

e^s = Be' [sinh (A a ss )] n , (52) 

The constant A does not appear explicitly in the final equations, but 
it enters into the calculation of A^ and 

During monotonic loading under steady state conditions the back 

• • 

stress and drag stress reach their saturated values, thus B and D are 
zero. Therefore, by equating equations (50) and (51) to zero and 
solving for the resulting steady state plastic strain rate gives 

l l ss = BqG 1 [sinh(A 1 B $s )] n , 


(53) 
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.1 

£ SS 


B C 0‘ 


sinh 


q SS~ B SS 

D SS 


1.5 ^ n 


, (54) 


.1 

E SS 


C 2 B c 0' [sinh(A 2 D ss 3 )] n 
C 2 + B SS " A 2 D SS / A 1 


(55) 


Thru separation of variables in equations (52) through (55), it may 
be seen that 


A 1 B SS 


°SS ‘ B SS 


n 1.5 


J SS 


A 2 D SS " Ao SS * 


Hence, from equation (56) explicit relations for the constants 
and A 2 may be obtained: 


and 


A a 

A 1 = ~B 


SS 


SS 



A 

1_B SS /a SS 


(57) 

(58) 


By linearizing equation (52) one obtains the relation 


ln(e^ s ) = n In [sinh(Ao^^)] + In (B^o') . (59) 

The constant A is chosen so that a plot of In (e^) versus 
In [sinh(Ao<-<-)] falls on a straight line. If the test data points 
fall on a straight line for a particular value of A, then Fig. 3 
shows that (CD/AB)/(EC/EB) = 1 , or 


In e 4 - In c 5 


ln[sinh(Aa 4 )] - ln[sinh(Aag)] 

In e-j - In 


ln[sinh(Aa-| )] - ln[sinh(Ao^)] 


. ( 60 ) 


Where , a^, and are the steady state asymptotic stress 
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obtained from hystresis loops executed at strain rates e^> 
and eg respectively. Once a value for A has been chosen, then 
from equation (59) the slope of lnU^) versus ln[sinh(Ao ss )] plot 
is n and the intercept is ln{B^ 0 ‘). 

From stress drop creep tests the value of B<- s can be 
established and the constants A-j and A^ may be computed from 
equations (57) and (58) respectively. In additon, once A^ and A 2 
are known, then = Dg may be found from equation (56). 

The constant in equation (50) governs the hardening rate 
of the back stress. This constant is analogous to the constants 
n.j + n 2 in Walker's model; therefore, the interested reader is 
referred back to Walker's model for a detailed explanation of 
their determination. 

The final two constants yet to be determined are H 2 and C,,. 

The parameter H 2 controls the mount of isotropic hardening 
produced by a given amount of strain. A best-fit value for 
can be tried and the resulting simulations compared with the 
test data for the same loading, to select the one which causes 
the closest agreement. The constant C 2 sets a "floor" under the 
drag stress and its value is estimated by a best-fit procedure to 
the actual test data. 

The material constants computed for Miller's model, simulating 
Hastelloy-X at 1800° F, were evaluated by Walker [42] and are 
tabulated in Table 4. 
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Table 4. Material Constants Used In Miller's 
Model For Hastelloy-X at 1800° F 


Miller's Constants 

Numerical Value 

E 

13.2 x 10 6 psi 

n 

2.363 

B c 0 ' 

2.616 x 10" 5 sec' 1 

H 1 

1 .0 x 10 6 psi 

A i 

1.4053 x 10" 3 psi' 1 

H 2 

100 psi sec 1y,n 

C 2 

5,000.0 psi 

A 2 

4.355 x 10' 12 psi' 3 

e'(t-O) 

0.0 

B(t=0) 

0.0 

D( t=0) 

8,642 psi 
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CHAPTER IV 

EVALUATION OF SELECTED THEORIES 

Evaluation of the unified constitutive theories considered 
herein is based upon numerical simulations of Hastelloy-X at 
1800° F in uniaxial form by specifying input strain histories 
and comparing output stress histories. Since this thesis addresses 
two separate problems, a section on the numerical time integration 
study is presented first, followed by a section detailing the 
material parameter sensitivity study. 

Numerical Time Integration Study 

Accurate integration of the constitutive relationship given 
by equations (1) through (4) is extremely important to the numerical 
solution of a nonlinear field problem. However, the integration 

■ r 

process '/is somewhat hindered because the coupled system of ordinary 
differential equations defining the state variables may be "locally" 
stiff and thus is sensitive to the time step size and numerical 
algorithm. The instability problems associated with integrating 
stiff O.D.E.'s can be minimized through proper selection of time 
step size, when using standard integration techniques, or through 
the use of specially tailored algorithms which check for stiffness 
and automatically take the appropriate action. 

In general, equation (5) may be integrated between time t and 


t + At by writing 
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e ! U+At) 


Ae = 


/ 


(t; 


£ = 


t+At 

/ e 1 dt 


(61) 


or 


A e 1 = e 1 (t+At) - e^t) = 


/ 


t+At 


e 1 dt , (62) 


where e 1 is defined by the particular constitutive theory being 
employed. 


Integration Algorithms Considered . Four of the more commonly 
used integration schemes were selected for use in this study. These 
include: Explicit Euler Forward Difference Method, Implicit 

Trapezoidal Method, Runge-Kutta 4th Order Method, and Trapezoidal 
Predictor-Corrector (Iterative) Method. 

The basic presumption of the Euler method is that the inelastic 
strain rate function is assumed to be a constant over an interval of 
interest. For all but the highly-nonlinear portion of the stress strain 
curve, this assumption is quite reasonable at constant strain rates. 
However, when the material response is nonlinear, the assumption may no 
longer be valid unless a sufficiently small time interval is chosen. In 
addition, the Euler method is sensitive to numerical round-off. The 
major advantage of using the Euler technique is its ease of implementa- 
tion. The approximation to the solution of equation (5) using Euler’s 
method may be written as 


A = At e ^ (t) 


(63) 



The Implicit Trapezoidal Method is another commonly used 
integration strategy from the family of single step schemes. 

Overall solution accuracy is improved by averaging the evaluated 
function over the desired interval assuming a small enough time 
increment has been selected. The Implicit Trapezoidal Method 
approximates the inelastic strain variable as 

Ac 1 = [ k l (t) + e 1 ( t+At) ] (64) 

The Runge-Kutta Method, which is also a single step scheme, 
provides a straight forward high order approximation for computing 
Ae 1 for various locations along the interval in question. The 
self-starting characteristic and small error per step of this 
algorithm has made it popular among researchers for numerically 
integrating differential equations. It is pointed out that, like 
most other numerical procedures for solving differential equations, 
the methods of Runge-Kutta are sensitive to variations in the 
selected interval. Too small a value leads to an excessive number 
of computations; on the other hand, too large of an interval may 
well lead to an answer that differs significantly from the true one. 
The 4th order Runge-Kutta integration method approximates the 
solution to equation (5) as 

Ae 1 = 1/6 [K-j + 2K 2 + 2K 3 + K 4 ] , (65) 


where 
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K 2 = At c 1 ( t + ijAt, e 1 ( t ) + ^K-j ) 

K 3 = At £ 1 ( t + %At, e 1 ( t) + hK 2 ) 

K 4 = At e*(t + At, e 1 ( t) + K 3 ) . 

The last integration algorithm selected for use in this study is 
the multi-step Trapezoidal Predictor-Corrector method. The objective 
of this strategy is to predict Ae* as accurately as possible using an 
iterative procedure. As a first guess, c*(t+At) is computed by the 
Euler method 


Predictor: e^(t+At) = e^tj + At c*(t, e*(t)). (66) 

This value is then substituted into the Trapezoidal scheme 

Corrector: c^t+At). = t) + y [e I (t,e I (t.) ) 

+ £*(t+At,e*(t+At) .)] •, (67) 

which may also be written as 


Ac j = Cc I (t) + e*(t+At)]. 


( 68 ) 


i. L 

The subscript j in equations (67) and (68) refers to the j 
iteration. If the values for Ae^ determined by (66) and (67) differ 
too greatly, equation (67) is re-evaluated using the most recent 
computed value of e*(t+At)j. This iterative process is continued 
until the convergence criterion, 


4e‘. 

J 


Ae^. ■, 

J +1 


- r tol * 


(69) 
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is satisfied. Selection of an appropriate r ^ is important 
because it governs solution accuracy, number of iterations per 
interval, and ultimate solution time. For this study r^ was 
chosen to be 0.000001. While the error of this iterative procedure 
is greatly reduced, it is still sensitive to numerical round-off. 

Strain Histories Considered . All numerical simulations were 
performed under strain control conditions using one of four input 
histories. These included: An alternating square-wave strain rate; 

a ramp up and ramp down strain rate; and two complex multi-step 
strain rate histories. Restrictions which were arbitrarily imposed 
on each loading sequence included: 1) the total response time 

would be 35 seconds; 2) the total strain would not exceed 0.5%; 3) 
all simulations would begin with a positive (or tensille) load; and 
4) the numerical experiments would end with a 10 to 15 second 
relaxation period. While these strain rate histories (see Fig. 4) 
were not developed to model any specific physical phenomena, they 
are representative of conditions under which a material could be 
strained either in the laboratory or the field. 

Qualitative and Quantitative Comparisons . The integration 
investigation was divided into two tasks. The first task included a 
comparison of the numerical stability of each constitutive theory 
in order to access the degree of nonlinearity of the prescribed 
growth laws. The comparisons were based upon predicted stress-strain 
behavior for constant strain rate tests simulating monotonic loading. 
Explicit Euler forward integration was used to obtain the stress 




1) Alternating Square Wave 
Strain Rate History. 



3) Complex Strain Rate 
History - 1 . 


Fig. 4. Input 


t (sec)' 



2 ) Ramp-up- ramp-down 
Strain Rate History. 



4) Complex Strain Rate 
History - 2. 



in Rate Histories. 
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histories presented in Fig.'s 5 through 8 for various time steps. 

From Fig. 8 it may be concluded that Miller's model is the most 
sensitive to step size variation because of the oscillatory response 
of the stress-strain curve for the larger At values. 

The oscillatory nature of Miller's model can be traced back 
to the hyperbolic sine function that is used to characterize the 
inelastic strain rate and other growth laws. In contrast, the power 
law function used by Bodner, Walker, and Krieg, to model a material's 
response appears to be much more stable for a wide range of time 
i ncrements . 

Figures 5 through 8 also reveal that as the step size is 
increased, all of the constitutive theories tend to overestimate 
the stress history in the 0.1 to 0.2 percent strain region and 
then automatically correct themselves as deformation continues. 

The "self-correcting" phenomenon of predicted stress appears to be 
an intrinsic property of the prescribed growth laws as is illustrated 
in Fig's. 9 and 10. Figure 9 is a plot of the inelastic strain 
computed by Walker's model for At's of 0.0005 and 0.1 seconds. If 
the smaller time step is defined to predict the correct results, then 
the inelastic strain computed when using the larger At is initially 
underestimated. Because of the coupled nature of equations (1) 
through (4), subsequent integration reverses this trend and results 
in an over prediction of the inelastic strain. While not shown in 
Fig. 9, the difference between the two curves assumes a more constant 
value as the total strain is increased. That is, the difference 
(or error) does not grow. Figure 10 shows the identical situation 
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Fig. 5. Monotonic Loading of Hastelloy-X at 1800° F 

Predicted by Bodner's Model (c = 0.001 sec -1 ). 
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Fig. 6. Monotonic Loading of Hastelloy-X at 1800° F 

Predicted by Walker’s Model (e = 0.001 sec' 1 ). 
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Fig. 7. Monotonic Loading of Hastelloy-X at 1800° F 
Predicted by Krieg's Model (c = 0.001 sec - l) 



Stress (KSI) 


ORIGINAL PAGE IS 
OF POOR QUALITY 


46 


At 

«* 

0 • S0CD-C3 

III 

cn 

; i. E R 

MODEL 

E L 1 L E R 1 

At 

MS 

0 -5000-02 

SEC 

MILLER 

MODEL 

1 EULER I 

At 

*= 

O -S00D-01 

SEC 

MILLER 

nODEL 

l EULER 1 

At 


0-1000 00 

SEC 

MILLER 

nOOEL 

1 EULER 1 

At 

S= 

0 • 1250 00 

SEC 

MILLER 

MODEL 

(EULER) 



Strain (in / in) 


Fig. 8. Monotonic Loading of Haste!loy-X at 1800° F 

Predicted by Miller's Model (e = 0.001 sec - l). 
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Percent Total Strain 

Fig. .9. Inelastic Stain Predicted by Walker's Model 
For Hastelloy-X at 1800° F. 
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arises in the prediction of the back stress. The over/under 
estimation of computed stress continues until the equations 
( (1) through (4) ) naturally dampen themselves out, as is shown 
in Fig. 6. It should be noted that the "self-correcting" nature 
of the constitutive models was only apparent when the Euler method 
was used to numerically integrate the prescribed growth laws. 

Figure 11 shows that when either the 4th order Runge-Kutta or 
Trapezoidal method was used, the peak values of stress were simply 
over- or under- predicted and there was no self-correcting tendency 
observed. 

The second task of the integration investigation was to 
evaluate the constitutive theories considered herein in terms of 
solution stability, accuracy, and computational efficiency when 
numerically integrated using various algorithms. A test matrix 
consisting of integration methods, time step sizes, and strain 
input histories was established in order to make qualitative and 
quantitative comparisons. However, before this phase of testing 
could begin both a baseline or "pseudo correct" solution for each 
input history was needed and the integration algorithms needed 
to be functionally verified so that, in the limit, they would 
predict the same response. 

For each strain history considered, a baseline or "pseudo correct" 
solution was obtained analytically, since no experimental data 
were available. The Euler method was used to numerically integrate 
the prescribed growth laws for each constitutive model until 
successive changes in step size predicted an identical response. 
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Fig. 11. Comparison of Integration Methods For Miller's Model For 
Equivalent Time Steps (Hastelloy-X at 1800° F). 
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This process was repeated for the various input histories, thus 
developing a family of baseline solutions. 

Verification of the integration algorithms was accomplished 
by using an a priori assumption that all numerical approximation 
procedures would yield the same output stress histories if the step 
size was small enough. Results of this investigation indicate 
that the various integration schemes converge to the same solution 
for the strain histories and constitutive theories considered 
herein. In addition, the computational efficiency of each method 
was obtained through selection of equivalent At's. Euler's method was 
observed to be the most efficient (in terms of execution time on 
an Amdahl 470V/8 computer), followed by the implicit Trapezoidal 
method (a factor of 2 slower), then the 4th order Runge-Kutta method 
(a factor of 4 slower), and finally the Trapezoidal Predictor- 
Corrector method (a factor of 12 slower). A plausible explanation 
for the relative slowness of the Predictor-Corrector method was the 
selection of an inappropriate convergence criterion. 

In the ensuing paragraphs qualitative and quantitative 
comparisons are made for the various test matrix parameters 
presented earlier. Since the amount of data generated during this 
study was voluminous, only a representative sample of the results 
will be presented. 

Figure 11 demonstrates the response predicted by Miller's 
model for a uniaxial, alternating square-wave strain rate input 
history. The growth laws were integrated using four different 
numerical approximations and a constant step size of 0.1 seconds. 
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The results indicate, for monotonic loading and unloading, that 
integration by the 4th order Runge-Kutta method consistently over- 
predicts the state of stress by as much as 15%. Conversely, the 
Trapezoidal method generally underestimates the stress history by 
4%. The Euler and Predictor-Corrector methods appear to provide 
accurate results at any point in the experiment. Similar observa- 
tions were made for Walker's, Krieg's, and Bodner's model, when 
subjected to the same strain rate history. These results indicate 
that for equal time steps, the 4th order Runge-Kutta method 
provides the least accurate results. 

Figures 12 and 13 illustrate that various constitutive models 
may differ appreciably when the same integration method is used 
(in this case the Euler method). In Fig. 12 it can be seen that 
Miller's theory gives a considerable oscillatory response for time 
steps greater than 0.05 seconds, whi 1 e Walker's theory shown in 
Fig. 13 gives a much smoother response for the same time step. 
Figures 12 through 15 also demonstrate that numerical stability 
and stress-overshoot are strain rate dependent. Note that in 
Fig's. 12 and 14 the oscillatory response of Miller's model does 
not occur until a strain rate of 0.001 in/in/sec is applied. 

In Fig's. 13 and 15 the amount of stress-overshoot (for the larger 
step sizes) appears to grow with increasing strain rate. 

Figures '16 and 17 present results predicted by Walker's model 
for a ramp-up-ramp-down strain rate history. Unlike the constant 
strain rate histories, this type of loading appears to affect the 
predictive capabilities of the model. As is illustrated, the 





Fig. 12. Stability and Accuracy of Euler Integration For 
Miller’s Model Using Various Step Sizes 
(Hastelloy-X at 1800° F). . 




Fig. 13. Stability and Accuracy of Euler Integration For 
Walker's Model Using Various Step Sizes 
(Hastelloy-X at 1800° F) . 
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tensile response is well behaved and the peak response is estimated 
fairly accurately for all step sizes. However, the compression 
response is oscillatory and inconsistent for At's greater than 0.05 
seconds. 

Figure 18 shows a 35 second response obtained with Bodner's 
theory using various integration algorithms. Time steps were 
selected to obtain solutions which required approximately equivalent 
computation times. As can be seen, the Euler method solution compares 
favorably with the baseline prediction, whereas the 4th order Runge- 
Kutta and Trapezoidal methods tend to over- or under- predict the 
peak response. Again, the time step used for the Euler method was 
smaller than those used in the other methods but the total computation 
time was the same. Similar observations were made for Walker's, 
Krieg's, and Miller's models when evaluated under the same conditions. 

Material Parameter Sensitivity Study 

Another important issue this research addressed was what amount 
of degradation in predicted response could be expected due to 
variations in material constants. Since each model has a number of 
material parameters which are determined through a set of complicated 
tests, the end user must be made aware of the effect small variations, 
i.e. typical experimental errors, have on the predictive capabilities 
of the constitutive theory in question. In addition, the experi- 
mentalist will better understand the importance of setting up and 
conducting accurate experiments from which the material constants are 


obtained. 
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Two procedures were developed to study the effect of material 
parameter variation. The first procedure characterizes the 
sensitivity of each constant within a particular 'constitutive 
model. Hence, the experimentalist will know what parameters need 
to be determined most accurately. Each constant’s sensitivity to 
variation was determined by adjusting it by a prescribed amount and 
comparing predicted results. 

The second procedure provided an "upper bound" on the effect 
that experimental error would have on the determination of the 
material parameters and ultimately, predicted response. Therefore, 
the experimentalist will know to what degree of accuracy the test 
data need to be obtained in order to insure accurate results. 

The effect of experimental uncertainty was studied by adjusting 
test data by a prescribed amount (simulating some type of error) 
and then using that data to determine the material constants. 

The overall material parameter sensitivity study was performed 
by subjecting Bodner’s, Walker's, Krieg's, and Miller's models to 
an alternating square-wave strain rate history. Euler's method was 
used to numerically integrate the prescribed growth laws using a 
time step of 0.005 seconds. 

All nine of the material constants in Bodner's model were 
individually adjusted by -5% in order to identify the parameters that, 
when changed, would have the most pronounced effect on the predicted 
response. Results of the sensitivity study indicate that the two 
most sensitive constants were the stress exponent, n, and the 
hardening recovery constant, D^, with the first of the two causing 
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the largest deviation from the baseline solution. The predicted 
response, after each parameter was changed by -5%, is shown in Fig. 19. 

The constant n, in Bodner's model, is determined from a series 
of constant strain rate tensile tests, where the stress has reached 
a saturated value. Through the use of equation (10) or by plotting 
the data as depicted in Fig. 1 (a), it can be seen that the only 
experimentally measured quantity that affects the value of n is the 
stress. The parameter D£ is obtained primarily from constant stress 
creep tests. Equation (13), which is used to compute D^, shows that 
the measured creep rate and the peak stress values (which are used to 
calculate the material constant n) can affect the determination of 
this parameter. Hence, the aforementioned experimental quantities 
must be measured very accurately in order to insure a minimal error 
in the predicted response. 

In addition to studying the sensitivity of individual constants 
in Bodner's model, the effect of a ± 5 % change in the test data 
(simulating experimental error) used to compute all of the material 
parameters was investigated. From the procedure discussed above, the 
stress exponent was shown to be the most sensitive constant. 

Therefore, the values of stress, determined from constant strain 
tests, were adjusted by ±5% in order to produce the largest change in n. 
Then the creep rates, from constant stress creep tests, were 
modified by ±5% in order to obtain the largest variance in D ^ • Hence, 
with new values of stress, creep rate, n, and D^, the other seven 
material parameters were re-evaluated and the resulting stress 
history, as shown in Fig. 20, was predicted. However, Fig. 20 
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indicated that as a result of the way the constants are determined 
in Bodner's model, a ±5% change in the test data used to compute 
the material parameters has little effect on the predicted response. 

The sensitivity of the nine constants used in Walker's model 
were evaluated by individually changing each parameter by -5% and 
comparing the predicted stress history to a baseline solution. 
Identified as the two most critical constants were the drag stress 
(assumed to be a constant value in this model), D, and the stress 
exponent, n, with the latter of the two being the most sensitive. 

Shown in Fig. 21 is the predicted response when each parameter 
(D and n) was adjusted by -5%. 

To compute both the stress exponent and the constant valued 
drag stress, the two maximum tensile stresses at the two highest 
strain rates (from steady state hysteresis loops) are needed. In 
addition to the above data, the saturated back stress in tension 
(B*; ), obtained from stress drop tests, is required. Values for 

n and D may then be determined from equations (30) and (31), 
respectively. In order to insure accurate output stress history 
predictions, precise experimental measurements of a and B^ ax are 
necessary. It should be noted that Walker [42] indicated that it is 
difficult to estimate the value of B^ ax because the creep rates are 
usually small and it is very difficult to discern when the creep 
rate ceases instantaneously. 

The effect that experimental error would have on the determination 
of the material parameters and ultimately, predicted response, was 
investigated for Walker's model. By adjusting the experimentally 
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measured values of stress from constant strain rate hysteresis loops 
and saturated back stress by ±5% to create a worst case scenario, 
each of the nine constants were sequentially recomputed. The effect 
of experimental uncertainty on the output stress history for Walker's 
theory is shown in Fig. 22. It can be seen that a ±5% error in all 
the test data requi red to compute the material parameters results in 
significant predicted response errors, with up to 30% over-prediction in 
the stress during the relaxation period. 

In a similar fashion. Fig. 23 illustrates the sensitivity of 

the stress exponent, n, and initial drag stress, D, to a -5% change 

in Krieg's model. Since these constants are analogous to the 

parameters, n and D, used in Walker's model, it is simply reiterated 

that the measured values of a and must be determined very 

max 

precisely in order to insure an accurate estimation of the predicted 
response. However, from Fig. 24 it can be seen that even after the 
constants are re-evaluated to reflect a ±5% change in the experi- 
mental test data, the output stress history is virtually unaffected. 
Hence, Krieg's theory appears to be more tolerant to experimental 
uncertainty than any of the models previously discussed. 

The nine constants in Miller's model were individually varied 
by -5% in order to characterize their sensitivity. However, unlike 
the preceding models , selectively changing any one parameter in 
this set of constitutive equations had little effect on the predicted 
response. In any event, the most sensitive constants were found to 
be n and A£. Figure 25 illustrates the effect that a -5% change in 
each of these parameters had on the computed stress history. 
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In order to minimize the error in the computation of n, the 
steady state asymptotic stresses obtained from hysteresis loops 
executed at several strain rates must be measured very accurately. 

In addition, to determine the parameter ^ one must precisely 
establish the value of the steady state back stress from stress 
drop creep tests. 

While Miller's model appears to be fairly insensitive to small 
variations in particular material constants, when the parameters 
were collectively re-evaluated (to simulate experimental error in 
the test data) an appreciable deviation from the baseline solution 
was noticed. Shown in Fig. 26 is the output stress history resulting 
from a ±5% change in the steady state stresses and back stresses 
that were used to recompute the respective constants. As can be 
seen, a significant over-prediction in peak stresses occurs. 
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CHAPTER V 
CONCLUSIONS 

In conclusion, the sensitivity of the unified constitutive 
theories proposed by Bodner, Walker, Krieg, and Miller, to numerical 
integration techniques and slight changes in material parameters was 
investigated. The integration sensitivity study involved the use of 
the Explicit Euler Forward Difference method, the Implicit 
Trapezoidal method, the Runge-Kutta 4th Order method, and the 
Trapezoidal Predictor-Corrector method. From the results of this 
study, the following observations were made: 

1) It appears that a simple integration scheme, like the Euler 
Forward Difference Method, is preferable to more complex 
schemes from the standpoint of accuracy, computation time, 
and ease of implementation. The results also indicate 
that under the conditions tested, the 4th Order Runge-Kutta 
method provided the least accurate results. 

2) The constitutive models may differ appreciably, in terms 

of numerical stability, when the same integration method and 
strain input history is used. Miller’s model was shown to 
be the most sensitive to time step variation, which is 
apparently due to the hyperbolic sine function that is used 
to characterize the inelastic strain rate function and other 
growth laws. 

3) The predicted stress response appears to be "self correcting" 
for large time increments when the constitutive equations 
are integrated by the Euler method. However, when either 
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the 4th Order Runge-Kutta or Trapezoidal method was used, the 
peak values of stress were simply over- or under- predicted 
and no self-correcting tendency was observed. 

The material parameter sensitivity study involved the 
characterization of the sensitivity of each constant within a 
particular constitutive theory and provides an "upper bound" on the 
effect that experimental error would have on the determination of the 
material parameters and ultimately, predicted response. From the 
results of this sensitivity study, the following observations were 
made: 

1) Small changes (5%) in the material constants may produce 
significant and unacceptable errors in the predicted stress, 
particularly the peak stress amplitudes during cyclic 
loading, and the amount of residual stress during stress 
relaxation. 

2) From. this numerical study, it would appear that more precise 
predictions can be achieved if the experimentally obtained 
data can be measured more accurately. The type of tests and 
data that are required to determine the material parameters 
for each model vary, however some generalizations can be 
made. For example, most theories require either monotonic 
or cyclic loading of the specimen in order to obtain peak 

or saturated values of stress for several constant strain 
rates. In addition, most models require a stress drop or 
constant stress creep test in order to determine the steady 
state or saturated back stress. While other specific 
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experiments may need to be performed, those mentioned above 
are probably the most critical to obtaining accurate 
analytical predictions. 

It should be noted that the conclusions reached herein are based 
on the use of uniaxial simulations. It is expected that the 
conclusions would hold for multi axial situations; however, further 
evaluation is needed in this area. 
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ABSTRACT 


The use of internal state variables in modeling of inelastic solids is 
gaining widespread usage in current research. Therefore, it is useful to con- 
struct a well-defined framework for internal state variable models which is 
based in continuum mechanics. The objective of this paper is to review and 
clarify the general theory of internal state variables and to apply it to in- 
elastic metals currently in use in high temperature environments. In this 
process, certain constraints and clarifications will be made regarding internal 
state variables. 

It will be shown that the Helmholtz free energy can be utilized to con- 
struct constitutive equations which are appropriate for metallic superalloys. 
Furthermore, internal state variables will be shown to represent locally av- 
eraged measures of dislocation arrangement, dislocation density, and inter- 
granular gracture. Finally, the internal state variable model will be demon- 
strated to be a suitable framework for comparison of several currently pro- 
posed models for metals and can therefore be used to exhibit history depen- 
dence, nonlinearity, and rate as well as temperature sensitivity. 


INTRODUCTION 

The prediction of inelastic behavior of structural materials at elevated 
temperature is a problem of great importance which has accordingly been given 
a great deal of interest by the research community in recent years. These 
materials exhibit substantial complexity in their thermomechanical constitu- 
tion. In fact, so complex is their material response that it could be argued 
that without useful a priori information, experimental characterization is fu- 
tile. The purpose of this paper is to show how the thermodynamics with in- 
ternal state variables can be utilized to emplace certain constraints on the 
allowable form of thermomechanical constitutive equations, thus providing some 
limited insight regarding experimental requirements. 

Historically, there have been two distinct approaches to the modelling 
of inelastic materials: 1) the functional theory [1], in which all dependent 

variables are assumed to depend on the entire history of independent variables; 
and 2) the internal state variable (ISV) approach [2], w’herein history depen- 
dence is postulated to appear implicitly in a set of internal state variables. 
Lubliner [3] has shown that in most circumstances ISV models can be considered 
to be special cases of functional models. For experimental as well as ana- 
lytic reasons numerous recently proposed models for the classes of materials 



discussed herein have been proposed in ISV form. Therefore, in this paper 
the ISV method will be reviewed as well as clarified and it will be shown that 
this general framework is useful in modeling metals at elevated temperature. 

The paper begins with a review of ISV theory, and this is supplemented 
with a section describing the procedure for constructing macroscopicaliv av- 
eraged internal state variables. These concepts are then applied to metals' 
at elevated temperatures. Finally, applications to boundary value problem 
solving techniques are discussed. 


REVIEW OF THE INTERNAL STATE VARIABLE (ISV) APPROACH 


The concept of internal state variables, sometimes called hidden vari- 
ables, was apparently first utilized in thermodynamics by Onsager [4,5] and 
numerous applications are recorded since the second world war [2,6-14]. Al- 
though not originally described for application to solids, the approach which 
will be discussed herein is due to Coleman and Gurtin [2]. 

In the theory of internal state variables applied to solids the follow- 
ing state variables are required in order to fully characterize the state of 
the body at all points x. and at all times t:* 


1) 

the 

displacement field 

“i * W 0 

; CD 

2) 

the 

stress tensor 

Ai - Vv c) 

; (2) 

3) 

the 

body force per unit mass 

f. - fyvo 

; (3) 

4) 

the 

internal energy per unit mass 

u = u(.X ,t) 

; (4) 

5) 

the 

heat supply per unit mass 

r = r(x k ,t) 

; (5) 

6) 

the 

entropy per unit mass 

s = s(x k ,t) 

; (6) 

7) 

the 

absolute temperature 

T = T(x k ,t) 

; (7) 

8) 

the 

heat flux vector 

q i = q i ( V C) 

; (8) 

9) 

k 

a. . 
ij 

= o£ (x ,t) , k = 1 , 2, . . 

ij m 

• , n 

; (9) 


where 0:^ • are a set of n internal state variables which are necessary to 
account for inelastic material behavior. Although they are listed here as 
second order tensors, they may be tensors of other rank as well [15]. 


* For convenience, only infinitesimal deformations will be considered here, 
although the general theory applies to finite deformations as well. 



The method of Coleman and Noll [16] may be used to obtain the spatial 
and time distribution of the body force f^ and heat supply r from the conser- 
vation of linear momentum and energy, respectively, assuming the displace- 
ments u^ and the temperature T are specified independent variables. Subse- 
quently, it is hypothesized that constitutive equations of state may be con- 
structed for the state variables described in (1) through (8) in terms of 
u. and T and their spatial derivatives: 

i 


i-J 


( X k>t) dj . (c. (x b , t ) , T(x^,t), g^Xj^jt), 


1J mn 


mn k 


u(x k’ C) = u(£ mn (x k jt) ’ T( V°, S m ( V C) ’ a mn (x k’ C)) 
s ( x k> t) = s(e mn (x k ,t), T(x k ,t), g m (x k ,t), aP n (x k ,t)) 

q i ( V C) = q i (£ mn (x k’ t} ’ T(x k’ C) ’ 8m ( V t} ’ a L (x k' t)) 


( 10 ) 

(ID 

(12) and 

(13) 


where g 

m 


is the spacial temperature gradient T, 


and 


ij 


%(u. . + u. ) 
i.J J. 1 


• (14) 


The form of equations (11) through (13) implies that all constitu- 
tive equations are evaluated in the specified state (x^.t). For this rea- 
son a-M , u, s, and q^ are termed observable state variables since they 
can be determined from equations of state even though there is implicit 
history dependence via the internal state variables , which are defined 
to be of the form: 



(£ , T, g , 

lj mn a m 


l x 

Ci ) 
mn 


where time and spacial dependence have been dropped for notational con- 
venience. If equations (15) are at all times integrable in time, then 
the following form is equivalent to (15): 


<.(x t) = 
ij m 


r 


D. . (x ,t’)dt’ 
ij m 


(13) 


(16) 


where t is the time of interest and t' is a dummy variable of integra- 
tion. Therefore, it is apparent that a^j are not directly observable 
at any time and must therefore be considered to be hidden or internal. 


Although the above framework has been shown to be applicable to rate 
dependent crystalline solids [17,18], it is often misconstrued that the ab- 
sence of explicit strain-rate dependence renders the model inappropriate 
for use in viscoplasticity theories. It is alternatively hypothesized that 


a. . 
ij 


0 . . (£ , 

ij mn 


£ , T, g , a ) 

mn m mn 


(17) 


is an appropriate form of thermomechanical constitutive equations (10). Al- 
though metals at elevated temperature certainly exhibit strain-rate depen- 
dence, there are several reasons why equations (17) are less desirable than 
equations (10). First, equations (17) are not actually equations of state 
since the inclusion of strain rate implies knowledge is required at some time 
other than the current time t. Secondly, as demonstrated in discussions of 



materials similar to (17) but without internal state variables [19], very 
little useful information will come from thermodynamic constraints. Finally, 
explicit strain rate dependence is actually redundant for the materials dis- 
cussed herein, as will be shown later. Therefore, although this is certainly 
a semantical issue, equations (10) through (13) and (15) are utilized as the 
constitutive model in the balance of this paper.. 


It should also be pointed out that internal state variable growth laws 
(15) could contain explicit strain-rate dependence: 


a. . = . ( e 

ij iJ 


ran* 


mn 


, T, g, a") 
mn mn 


, (18) 


as in the example of a rate independent elastic-plastic material , in which 
equations of the above form are linear in strain rate: 


• Z n k , 

a. . = . (e , 

ij ljpq mn 1 


T, 


A .• 
a ) e 
mn pq 


(19) 


Such a form, although not excluded by the principle of equipresence [20], 
is only necessary in the circumstance wherein specific rate independence is 
required, as can be demonstrated by direct substitution of (19) into (16). 
Furthermore, although the thermodynamic constraints will vary somewhat when 
(19) are utilized [21,22], the results will be quite similar to those described 
below. 


On the basis of the Coleman-Mizel procedure [23] it can be shown that 
satisfaction of the first and second laws of thermodynamics for the class of 
materials detailed above will lead to the following conclusions: 

h = u - Ts = h ( e , T, c£) ; (20) 

mn mn 

where h is the specific Helmholtz free energy; 


'k£ 


s 



3h ‘ 

3T 


; ( 21 ) 
; ( 22 ) 


and 


q i = _k ij 8 j + ° (8 i ) ’ (23) 

Equations (21) should not be interpreted as defining as hvperelastic material 
since the Helmholtz free energy, described by (20), is dependent on the in- 
ternal state and therefore path dependent. 


Although not directly related to our problem, it is useful to note that 
the path dependence of the Helmholtz free energy precludes the usefulness 
of equations (21) in Rice's J-integral for fracture mechanics [24]. However, 
in the case wherein the loading path is radial: 




i = vr - ~ 

ij ij 


, ( 24 ) 



where kjj and kj_j are constant coefficients, then it is well known that equa- 
tions (15) are directly integrable so that the free energy can be described 
by 


h * h(£ 


tan 


T, a ) = h(e . 
pq mn 


T, a (£ )) = h(£ T) 

pq mn mn 


• (25) 


Thus, for the case of proportional loading only, the constitutive equa- 
tions are derivable directly from a potential function and the J-integral 
method is applicable. 


THE LOCAL AVERAGING PROCESS 


Constitutive equations (10) through (13) and (15) are theoretically 
pointwise in nature; that is, they are applicable to fixed infinitesimal 
material points. However, practically speaking, there is no way to construct 
experiments on material points since at the microscopic level the continuum 
assumption becomes invalid. Rather, it is considered acceptable to construct 
constitutive equations by subjecting local specimens to surface deformations 
(or tractions) which lead to spacially homogeneous stresses and strains so 
that some local average of the pointwise observable state variables can be 
determined directly from the effects on the boundaries of the specimens. 


As shown in Fig. 1, the scale of the smallest dimension of a local speci- 
men is generally constructed so as to be at least an order of magnitude larger 
than the scale of the largest material inhomogeneity. This sizing helps pre- 
serve the continuum assumption while at the same time averaging out the effects 
of point defects such as crystal lattice dislocations. Conversely, the scale 
of the largest dimension of a typical specimen should be as small as possible 
compared to the scale of the global boundary value problem of interest. This 
constraint is necessary in order to pre- 
serve the notion that constitutive equa- 
tions are indeed pointwise in nature, but v __ 

it is pragmatic in that it is a simple 

matter of economy. □ 


The local rather than pointwise con- 
stitutive equations that result from ex- 
perimentation are assumed to be of the 
same form as pointwise equations (10) 
through (13) and (15). For example, in 
the uniaxial test described in Fig. 1 
it is customary to define 


global 
domain of 
interest 


Y 


‘! 


■i \ 


i 


u 


'u 


_i 

A 


/ 

B 1 

/■ 


°1 i dx 2 dx 3 


(26) 


physical 
inhomogeneitied 


tzx 


local 

specimen 


L 

i 

± 


(27) 


Figure 1 



and 


, ( 28 ) 


T - T (a ^ , a^> a^) 

where L is the local specimen gage length, A is the cross-sectional area in 
the X9-X3 plane, and (ap a 2 , a 2 ) is some arbitrary point on the surface of 
the specimen. Utilizing these quantities, it is then hypothesized that 

o ,(£,,, T, ct\) = 0 ( e , T, af ) , 0 

11 11 mn 1111 mn 


where 


-i _ 1 / l . . 

a = 77 I a dx. dx„ dx,, 

mn V J mn 1 2 3 


and all quantities with bars represent the locally measured state variables. 

Although equations ( 29 ) represent an often used way of relating point- 
wise equations to experimental results, the local averaging process is never- 
theless fraught with shortcomings since definitions ( 26 ) through ( 28 ) all re- 
resent nonunique relations between pointwise state variables ^ij , £ij > I> 
and their locally defined counterparts O^j , , T, and . There are in 

fact an infinite number_of distributions a^(xj, x 2> x 3^ which will result 
in identical values of a^. However, assuming that the scale of inhomogen- 
eities is small and that the distribution of 0^ is random the specimen will 
be statistically homogeneous and the relation between a~ n and a^ n will be 
reasonably one to one. 


For example, suppose that during some monotonically ^Increasing local 
strain history a particular internal state variable cq 1 such as a single 
dislocation arrangement is governed on a pointwise basis by the almost dis- 
continuous behavior shown in Fig. 2 . Suppose further that the time t at 
which the internal state begins to change 7 
is determined by the pointwise stress Kf 

state. Then the number of dislocation 
rearrangements occurring in the local 
specimen as a function of time might be 
distributed as shown in Fig. 3 . If the 
local specimen is large compared to the 
scale of the dislocation, and there are 
numerous dislocation rearrangements, as 



is usually the case in testing of metals, 
then the peak of the curve shown in Fig. 

3 will be several orders of magnitude 
greater than unity. It follows from 
equations ( 30 ) that the locally averaged 
value of the internal state variable 
represented in Fig. 2 will be as quali- 
tatively shown in Fig. 4 . 



Figure 2 



APPLICATION TO METAL CONSTITUTION 


In order to describe the class of metals discussed herein, the free en- 

p 

ergv maybe expanded in terms of the elastic strain tensor and the temper- 
ature T in a second order Taylor series expansion as follows: 



, (31) 


T I 

where is the thermal strain tensor and is the inelastic strain tensor, 
considered to be an internal state variable [17,18,22,25], and 


h 


h R + 2p 


F p C 

„E E _J£ 

°ij ijklScl “ 2T 


(T - T r ) 


, (32) 


where the subscript R refers to quantities in the reference state, ^ is 


is the specific heat at 


the linear elastic modulus tensor, and C v 

constant elastic strain. Substitution of equation (32) into (21) will result 
in 


= _ T /l i 


O. . = D, i (£ 
kl klmn mn 


I T 

s - e 1 ) 

mn mn 


(33) 


The above equations, together with internal state variable growth laws (15), 
will be shown to be a suitable framework for comparison of all of the models 
to be discussed herein. 


Internal State Variables in Metals 


It is now generally agreed in the literature that in single crystals there 

are two locally averages internal state variables: the back stress (a ) 

lij 

representing dislocation arrangement; and the drag stress (&->) representing 
dislocation density; where the bars have been dropped for convenience and the 
superscript has been converted to a subscript in order to avoid the confusion 
which would arise if a state variable were raised to some power. For obvious 
reasons the back stress is a second order tensor, whereas the drag stress is 
a scalar. In specimens composed of multiple crystals it is generally agreed 
that a third internal state variable loosely termed damage (a^ ) is neces- 

ij 

sary in order to account for intergranular mechanisms such as grain boundary 
sliding and microvoid growth and coalescence that may occur at high tempera- 
ture and/or large strain. Although damage is obviously a directionally re- 
lated quantity and therefore tensorial in nature, it is difficult to distin- 
guish phenomenologically between damage and drag stress since both are pri- 
marily stiffness reducing mechanisms. 

Within the thermodynamic framework described earlier it is also possible 
to define the inelastic strain tensor to be an internal state variable. How- 
ever, this interpretation is not generally utilized within the materials 



literature. It is hypothesized that the 
rate of growth of the internal state vari- 
ables does not depend on the inelastic 
strain tensor so that 


no. of dislocation 
rearrangements in 
a local 4 
specimen 


*ij = Q ij (£ mn’ T ’V a i > a 2’ a 3 } 

mn mn 


(34) 


Due to the form of equations (34) it is 
said that since the inelastic strain ten- 
sor does not appear on the right hand 
side it is not an internal state vari- 
able. However, within the framework de- 
fined herein, it is still possible to 
construct an internal state variable 
growth law of the form 



Figure 3 


*lj E n ij (£ mn' T ’V a l 'V“3 } 

mn mn 


(3; 


which is precisely in agreement with def- 
initions (15) . 
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Figure 4 


In order to qualitatively verify the 
supposition that the inelastic strain ten- 
sor can be regarded to be an ISV, consider 
the example of a uniaxial bar subjected to 
applied displacements such that the end 
tractions will be evenly distributed. It 
is customary to deduce the inelastic strain 
in an experiment of this type by utilizing 
the output from a load cell to determine 
the stress and then making use of equations 
(33) to determine the elastic strain. This result and the total strain mea- 
sured by an extensometer are then substituted into equations (31) to deter- 
mine the inelastic strain. Nevertheless, this does not imply that the inelas- 
tic strain tensor in an observable state variable. This result can be arrived 
at only in constitutive experiments such as uniaxial bar tests in which the 
stress and strain tensors are spacially homogeneous. In heterogeneous boun- 
dary value problems, only two state variables may be input (temperature and 
either stress or strain), and for this case equations (31) and (33) must be 
supplemented with an ISV growth law of the form of equations (35) in order 
to determine the inelastic strain tensor. Therefore, in the context of the 
current thermodynamic framework the inelastic strain tensor may be interpreted 
to be an ISV. 


A Framework for Current Metals Models 

In order to establish that current models can be constructed from equa- 
tions (33), consider the standard solid shown in Fig. 5. The governing dif- 
ferential equation for this analog is 

^1 • E co 

o+ — o = E„£ + n M [i + 


+ p— 0 = E 

e m 




, ( 36 ) 



where by convention the stress is denoted 
G and the strain is £. Equation (36) may 
be written in the following equivalent 
form 


!V E ~> 




(37) 


In accordance with the instantaneous lin- 
ear elastic behavior of metals, it is as- 
sumed that 


E+E m - E = Young's modulus = constant, (38) 

so that it is clear that equation (37) can 
be integrated in time to give the following 
stress formulation 


a(t ) 

e(t i ) = ~r~ + e (t i } • (39 > 


where £ is the inelastic strain, defined 
by 


e I (t 1 ) 



(40) 



f 


f 

a 


Figure 5 


Equation (39) may be solved for the stress and substituted into equation (40) 
so that it is clear that equation (40) is in accordance with ISV growth laws 
(16). Further, it can be seen from the standard solid analog in Fig. 5 that 
since O-E^e represents the stress in the Maxwell element, is not observ- 
able, so that satisfies the two conditions required for it to be an internal 
state variable. 


Equation (39) may be written equivalently in the following strain for- 
mulation: 

a(t L ) = EUUp - e I (t 1 )] , (41) 

which is an equation of state compatible with constitutive equations (10) as 
well as equations (33) . Since no other internal state variables are present 
in this equation, and also, no additional internal state variables are present 
in growth law (40) it is apparent that the standard solid analog with constant 
coefficients E^, and E^ is a single internal state variable model. 

It has been noted by several researchers that the standard solid is an 
appropriate analog for thermoviscoplastic metals if the springs and dashpot 
are nonlinearized [26,27]. In order to demonstrate this feature, consider a 
multiaxial extension of equation (36): 



, ( 42 ) 


O+K o = G £ + M £ 

pq pqmn ran pqmn mn pqmn mn 

where by convention the small strain tensor £^j is used in conjunction with 
the work conjugate stress tensor G-jj . In order to model metals K-pqjnn , G pqmn , 
and Mpq mn are required to be nonlinear in some as yet undetermined way. In 
addition, in accordance with constraint equation (38), it is required that 


K. 1 M . = D. . 

ljmn mnkl ljkl 


> (43) 


where D^j^l is the linear elastic modulus tensor. Equations (42) may be re- 
written in a strain formulation equation of state form as follows: 


o . . 

ij 


D [> -si 

ijkl L 'kl kl 


where 


J. 

°kl 


I 

£ . . 


is the inelastic strain tensor, defined by 

M.* [a -G £ ]dt 
ljpq pq pqmn mn 

_co 



(44) 


(45) 


Substituting equations (43) and (44) into equations (45) will result in 




(K . ^ [ £ 

ljmn mn 


z l ] - mT l 

mn ljpq 


G £ } dt 

pqmn mn 


, (46) 


so that equations (46) are in accordance with growth laws (16). The number 
of internal state variables contained in the model will depend on the degree 
of nonlinearity proposed in the nonlinear tensors Kp qmn , G pqmn , and M pqmn , 
and this will be discussed in the following section. However, before continu- 
ing, it should be pointed out that the constitutive equations developed in 
this section assume that the elastic and inelastic strain tensors may be lin- 
early decoupled. It has been shown that this assumption is invalid for finite 
deformation [28]. However, even under finite deformation conditions the in- 
elastic strain is decoupled from the elastic strain in such a way that the 
inelastic strain tensor may be considered to be an internal state variable. 


Current Models for Metals 


The framework for metals models discussed in the previous section can 
be used to describe numerous models currently under development [26,27,29-58]. 
For example, the microphysicallv based isothermal model proposed by Krieg, 
et al. , [30] is of the form described by equations (33): 


= D ( £ - p 1 ) 

ijkl v kl kl' 


( 47 ) 


o . . 
ij 


where 



, ( 48 ) 



and £ 0 and m are material constants, and is the deviatoric stress tensor 
and a[^. is the deviatoric component of the back stress tensor. Since equa- 


tions (48) contain the stress tensor, substituting equations (33) into (48) 
will result in equations consistent eith growth laws (15). In addition, Krieg, 
et_al. , give the back stress and drag stress to be, respectively, 


a 


a = A £ . . 
1 . . aii 

ij 


1 . . 


a ( a f a i f 

\ pq pq/ 


, (49) 


and 


“2 -*R < 


•I h 
£ 7 . £ 7 ./ — r D 

ij i] R 


, (50) 


where A a and A^ are hardening constants, and r a and r R are recovery functions 
of temperature and internal state variables. It can be seen that since ISV 
growth laws (49) and (50) are consistent with equations (15) , the model pro- 
posed by Krieg, et al . , contains three internal state variables: the inelastic 

strain tensor, the back stress tensor, and the drag stress tensor. 


Furthermore, classical plasticity theories can be described by the gen- 
eral form 



D, . (£ - E 1 ) 

ijmn mn mn 


where 


, (51) 




• (52) 


A is a scalar valued function of state, and F is a scalar valued state func- 
tion for inelastic behaviour often taken to be the yield function. If F is 
described by the von Mises yield criterion [53], given by 


F(a ii " a i > = ^< a ii ~ a i ) ( a ii - ) = k 2 , (53) 

ij ij J ij 

where is a tensor describing the yield surface center in stress space and 

ij 

k is a constant representing the yield surface size, then equations (52) can 
be written as 



, ( 54 ) 


e. . 

ij 




) 


resulting in a kinematic hardening model with constant yield surface size. 
Substitution of equations (51) into the above will yield a result consistent 
with rate independent ISV growth laws (19). 


Furthermore, if the yield surface translation is derived from the Ziegler 
modification [60] of the Prager work hardening rule [61], it may be described 
by 


“l = ^ (a ii " a i * ’ ( 55 ) 

ij J ij 

where p is a scalar valued function of state. By use of equations (51), equa- 
tions (55) can also be shown to be consistent with equations (19). Therefore, 
a classical plasticity-based kinematic hardening model contains two internal 
state variables: the inelastic strain tensor and the yield surface transla- 

tion tensor representing the back stress. 

In order to further illustrate the applicability of equations (33), (35) 
and (15) to current models for metals, ten of these models have been cast in 
uniaxial form in Table 1, wherein it is shown that although the framework for 
each model is identical (Valanis 1 model is in simplified form), the ISV growth 
laws vary widely both in number and form. 


CONCLUSION 


The main content of this paper has been to review and clarify the continuum 
and thermodynamics based internal state variable model for application to ther- 
moviscoplastic metals. In this process the following points have been made: 

1) the definition of an internal state variable utilized in this model 
has been clarified; 

2) internal state variables in metals represent local averages of dis- 
location arrangement, dislocation density, and intergranular damage, 

3) in the context of the ISV definition given here, inelastic strain may 
also be interpreted as an internal state variable; 

4) the path dependent Helmholtz free energy may be expanded in a second 
order expansion in elastic strain and temperature in order to obtain a stress- 
strain equation of state; 

5) rate dependence enters the constitutive equations implicitly via the 
inelastic strain, as demonstrated by the nonlinear standard solid analog; and 

6) a three-dimensional generalization of the standard solid may be used 
as a means of comparison of the general form of several currently proposed 
models . 



Further ramifications of the ISV model discussed are also of importance, 
although not detailed herein. For example, this model may be utilized to con- 
struct a coupled heat conduction equation which may be utilized to predict 
heat generation in thermoviscoplastic metals [62]. Furthermore, the concept 
of internal state variables may be utilized to construct models for the mech- 
anical constitution of composites with damage [63,15,65,66]. 
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ABSTRACT 


A Comparison of Current Models for Nonlinear Rate-Dependent 
Material Behavior of Crystalline Solids. (May 1986) 
Joachim Michael Beek, B.S., Texas A&M University 
Chair of Advisory Committee: Dr. D. H. Allen 

This thesis reviews three theories for prediction of 
inelastic deformation which are based on considerations of 
the microstructural behavior of materials. These theories 
are those of Krieg, et al . , Bodner, et al . , and Miller, 
et al . The thesis opens with a review of the mechanics of 
continua with internal state variables and a review of the 
historical development of constitutive modelling. A 
detailed discussion, including chronological development, 
review of theory, and method of determination of material 
parameters, of each model follows. An experimental data 
base. is established from which the material parameters of 
the constitutive equations are calculated. Finally, 
computer simulations of various load histories are 
performed and compared to experiment, and conclusions of 
the ability of each theory to model inelastic deformation 


are drawn . 
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INTRODUCTION 


The prediction of inelastic material behavior in metals 
is a problem of great importance which has accordingly been 
given a great deal of attention by the research community 
in recent years. Recent technological advances require 
materials to function in severe chemical, mechanical, and 
thermal environments such as nuclear reactors and ga3 
turbines. Superalloys have been developed that are able to 
sustain loads at extreme temperatures; however, due to this 
environment, these materials exhibit substantial complexity 
in their material constitution in that they are highly 
nonlinear, rate-dependent, temperature-dependent, and 
history-dependent. 

Numerous theories have been proposed to predict the 
thermomechanical behavior of inelastic solids at elevated 
temperatures. These theories are usually based on one of 
the following concepts: 1) rate-dependent extensions of 
classical plasticity theory retaining the concept of a 
yield surface, 2) thermodynamics and/or nonlinear 
viscoelasticity theory, and 3) considerations of the 


The style and format of this thesis follow the Journal 
of Engineering Materials and Technology . 
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microstr uctural behavior of the material. It is important 
that each of these theories be able to model diverse 
phenomena associated with inelastic deformation including 
anelasticity (completely recoverable time-dependent strain 
at zero load), the Bauschinger effect, cyclic strain 
harden ing /sof ten i ng , rate sensitivity, creep, relaxation, 
and tensile and compressive loading and unloading. These 
are illustrated in Figures 1 and 2. 

The purpose of this research is to review three 
theories based on considerations of microstructural 
behavior. Since each theory is proposed for a material 
and temperature of its authors’ choosing, a set of 
experiments is performed by .this author with a single 
candidate material in order to create a common experimental 
data base with which the mathematical representation of 
each model can be characterized. The theories are 
subsequently analyzed qualitatively and quantitatively over 
a range of material behavior to determine their suitability 
in predicting inelastic deformation. 

This thesis opens with a short review of the mechanics 
of continua with internal state variables. This is 
followed by a literature review and, utilizing the 
framework of internal state variables, a detailed study of 
the models chosen for discussion. The experimental program 
associated with the research for this thesis is presented 
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as an assessment of the experimental data base requirements 
for each model, accompanied by the results of the 
experiments. A uniaxial comparison of the predictive 
capabilities of the models follows and, finally, 
conclusions from the complete survey are presented. 



Creep 


Anelasticity 


Fig. 1 Phenomena associated with inelastic deformation 
(uniaxial loading) 
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THE EQUATION OF STATE APPROACH 


The concept of internal state variables, sometimes 
called - hidden variables, is gaining widespread usage in 
current research on the modelling of inelastic solids. 
Therefore, it is useful to construct a well-defined 
framework for internal state variable models which is based 
on continuum mechanics. 

The concept of internal state variables was apparently 
first utilized in Onsager's work in thermodynamics [1,2] 
and has found numerous applications since the Second World 
War [3-12]. 

In the theory of internal state variables applied to 
solids, Coleman and Curtin [3] defined the elastic field 
problem to be one in which the following state variables 
are required in order to characterize fully the state of a 
body at all points and at all times t: * 

1) the displacement field u^ = u^. (x^.t); (1) 


* For convenience, only infinitesimal deformations will 
be considered here, although the general theory applies to 
finite deformations as well. 
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2) the stress tensor 


°ij " °ij (x k ,t); (2) 


3) the body force per unit mass = f^(x ,t); (3) 


) the internal energy per unit mass 

u = u ( Xj , t ) ; ( i|) 


5) the heat supply per unit mass r = r(x.,t); (5) 

J 


6) the entropy per unit mass s = s(x.,t); (6) 

J 


7) the absolute temperature T = T(x.,t); (7) 

J 


and 


8) the heat flux vector 


q j[ = q jL ( x^. , t) . (8) 


Thermodynamic constraints will lead to the conclusion that, 

for elastic bodies, entropy is generated only through heat 

conduction [33. For inelastic bodies, however, this 

conclusion cannot be drawn because additional entropy is 

generated. The state cannot be character i zed solely from 

the observable state variables; rather, a set of internal 

state variables a. . is also necessary, and these are 

^ J 

determined by observation of the entire past history of 
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states of the body: 

a l j = a i j ( x m’ t) ’ (9 ) 

where the superscript k ranges from one to the number of 
internal state variables necessary to characterize fully 
the state of the body. Although listed here as second 
order tensors, they may be tensors of other rank as well 

[ 133 . 

The method of Coleman and Noll [1-43 may be used to 
obtain the spatial and time distribution of the body force 
f £ and the heat supply r from the conservation of linear 
momentum and the conservation of energy, respectively, 
assuming the displacements u^ and the temperature T are 
specified independent variables. It is hypothesized 
subsequently that constitutive equations may be constructed 
for the remaining state variables in terms of u . and T and 
their spatial derivatives: 


°ij (x k ,t) ~ 


°i j ( e mn ( x k ’ ,T(x k’ ’ S m ( x k* t} * a mn ( x k ’ t} } ; (1 0) 
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u(x k , t) = 


U<£ n,n <x k' t) - T( V tJ 


8 ,» ( V t, ’“mn ( V U) 


(11 ) 


s( x k , t) = 


sU mn < V t) - T(x k 


t) 


8 n, <x k' t) '“,iin < V t>) 


( 12 ) 


q i (x K' t) 


q i (c mn (x k- t) ' T(x k' t) * g m (x k’ t) ‘ a mn (x k' t)); (13) 


where g m is the spatial temperature gradient T , m and 


e = - ( u . ,+u. .) 

■ij 2 i,j j,i 


(1 4 ) 


is the infinitesimal strain tensor. The form of equations 

(10) through ( 13 ) implies that all constitutive equations 

are evaluated in the specified state x k and t. For this 

reason, o , u, s, and q. are termed observable state 
^ J i 

variables since they can be determined from equations of 

state, even though there is an implied history dependence 

through the internal state variables a* 3 . These are 

mn 

defined to be of the form: 



0 5) 


1 0 


ij 


0 . . ( e , T 
ij mn’ 


m 


«Sn> 


where time and spatial dependence have been dropped for 
rotational convenience and the superscripts k and p range 
from one to the number of internal state variables 
necessary to characterize fully the state of the body. If 
equations (15) are integrable in t at all times, then they 
can be rewritten as 


a* (x.t) = fljl k A x , t ’ )dt ' , (16) 

i j m •» ij m 

where t is the time of interest and t' is a dummy variable 
of integration. 

The above framework has been shown to be applicable to 
crystalline solids [15,16] and further discussion of this 
subject is provided in reference [17]. 

On the basis of the Coleman-Mizel procedure [18] it can 
be shown that satisfaction of the first and second laws of 
thermodynamics will lead to the following conclusions: 

h - u - TS - HU mn ,T,„P n ), (17) 

where h is the Helmholtz free energy; 


°ki = p dh/6e ki : 


( 18 ) 



1 1 


s = -dh/dT; 


( 19 ) 


and 


q t - -K ljgj * 0(g). 


( 20 ) 


where k is the thermal conductivity tensor. Although 
further results are obtainable [3. 15-17], they are not 
pertinent to the current research. The importance of the 
results above is that, in order to construct a complete 
description of constitutive equations (10) through (13) and 
(15), it is necessary only to prescribe the Helmholtz free 
energy. In order to describe the class of materials 
discussed herein, the free energy is expanded in terms of 

E 

the elastic strain tensor and the temperature T in a 

second order Taylor series as follows: 


h = h + e E .D. .. . ef , /2p 
R ijijklkl K 


c y (T - T r )/2T, 


(21 ) 


where the subscript R refers to quantities in the reference 

state, D. .. . is the linear elastic modulus tensor, 
ljkl 


V 


- T ( d 2 h/dT 2 ) is the specific heat at constant elastic 


volume, and 


C E I T 

= e . . ~ £ . . - e.. 

1J -tj ij 


( 22 ) 







The above equations, together with internal state variable 
growth laws (15), will be shown to be a suitable framework 
for comparison of all the models discussed herein. 

Theoretically, constitutive equations ( 10 ) through (13) 
and ( 15 ) are applicable to fixed infinitesimal material 
points. In practical terms however, it is not possible to 
construct experiments on material points since the 
continuum assumption becomes invalid at the microscopic 
level. Rather, it is considered acceptable to construct 
constitutive equations by subjecting local specimens to 
surface deformations (or tractions) which lead to spatially 
homogeneous stresses and strains. In this manner, a local 
average of the pointwise observable state variables can be 
determined directly from the effects on the boundaries of 
the specimen. 

As shown in Figure 3, the scale of the smallest 
dimension of a local specimen is generally assumed to be at 
Least one order of magnitude larger than the scale of the 





largest material inhomogeneity, while the scale of the 
largest dimension of the local specimen is much smaller 
than the scale of the domain of interest. The former 
preserves the continuum assumption, while the latter 
conserves the notion that constitutive equations are indeed 

i 

) 

pointwise in nature. The local constitutive equations that 
result from experimentation are assumed to be of the same 
form as pointwise equations (10) through ( 13 ) and (15); for 
example , 

•„ ' «n<'n- T -“Sn>- (21|) 

where 

”11 • i ; B, , 11 (1 V , 3' <25) 

E n " l^l e i i dx i ' ' * 26 ; 

and 

T = T ( a 1 ,a 2 ,a ) , (27) 

where L is the length of the local specimen, A is the 
cross-sectional area in the x^-x^ plane, is an 


arbitrary point on the surface of the specimen, and 



( 28 ) 


1 5 


>< p « f. /„a P dx,dx_dx_, 

mn V V ran 1 23 


It Is then assumed that equation (2U) is a suitable 

replacement for equation (10). There are obvious 

shortcomings in this assumption in that equation (2^0 

represents a solution to a boundary value problem. This 

implies volume averaging of local balance laws, which 

violates the local nature of constitutive equations. In 

addition, there may be a question of uniqueness resulting 

from the fact that more than one global state of a given 

internal state variable will map into the globally averaged 

state variable described by equations (28). However, 

assuming that the scale of inhomogeneities is small and 

that the distribution of a P is random, the specimen will 

mn 

be statistically homogeneous and the relation between a P n 
and a P n reasonably one-to-one. The above procedure 
remains the only reasonable method for constructing 
constitutive equations from experimental data bases. 
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LITERATURE REVIEW AND GENERAL DISCUSSION OF MODELS 


Man’s curiosity about his environment and his constant 
striving to harness the power and resources of nature with 
increasing efficiency have been the driving force behind 
the study of the laws of nature and the attempts to use 
these laws to predict the solutions of many problems. 

The theory of elasticity dates to the 17th century with 
Robert Hooke's discovery that the tension in a spring is 
directly prop or t ional to its extension. This heralded the 
beginning of the characterization of mechanical behavior of 
materials using mathematical constitutive models . 
Contributions to the theory of elasticity were also made by 
Galileo, the Bernoullis, Navier, Cauchy, Euler, and 
Sa int- Venant , to name a few [22]. 

Classical Plasticity Models 
Often only a small part of material deformation is 
elastic and, in general, a material will not behave 
elastically at large strains. Consequently, in 1864 Tresca 
[23] proposed a criterion to predict the onset of inelastic 
deformation in metals based on the maximum shear stress 
reaching a critical value. This was followed in 1913 by von 
Mises [24] with a yield criterion based on an equivalent 



stress state reaching a critical value; this has been 
interpreted as the elastic shear strain energy reaching a 
critical value [25]. 

Independently of each other. Levy [26] and von Mises 
[24] formed the basis of classical plasticity by proposing 
a three-dimensional relationship between the strain 
increment and the stress tensor which became known as the 
Levy-Mises equations: 


de 


ij 


d A 


ij 


(29) 


where de.. is the total strain increment, dA is a scalar, 
ij 

and ojj is the deviatoric stress tensor, given by 



°ij 


3°kk 6 ij 


(30) 


where 6^ is the Kronecker delta. In equations (29), the 
elastic strain component is assumed to be negligible; these 
equations consequently describe a fictitious, rigid, 
perfectly plastic material. Equations (29) were extended 
to take into account the elastic strain component by 
Prandtl [27] in 1924, who solved the plane stress case, and 
by Reuss [28] in 1930, who generalized the equations to 


three dimensions: 







where de^ is the plastic strain component and dX is a 
scalar history-dependent material property. The total 
strain increment de^ is then assumed to be the sum of the 

P 

plastic strain component de^ and the elastic strain 
£ 

component de ^ j . 

These developments enabled the classical plasticity 
theory to model elastic and rate- independent inelastic 
deformation fairly accurately for most metals with only a 


few restrictions. These were that the material be 


subjected to monotonically increasing loads and low 
non- transient homologous temperatures; that is, 
temperatures less than three tenths of the absolute melting 
temperature. At- higher temperatures, however, rate 
dependence of material behavior becomes significant, 
limiting the use of the existing plasticity theory. 

Rate dependence was initially treated as an extension 
of classical plasticity theory by Bingham [29] in 1922 by 
using the concept of a yield function dependent on the 
excess of the stress intensity over the yield stress for 
the case of simple shear. This was generalized by 
Hohenemser and Prager [30] in 1932 for the 
three-dimensional case. Further progress in plasticity 
theories was slow and advanced vastly and diversely only 


after completion of World War II. Prager [31] and Ziegler 
[ 32 ] proposed work-hardening rules for rate- independent 
plasticity to allow the yield surface to translate to model 
kinematic hardening. Drucker ' s stability postulate [333 
serves as a basis for the incremental theories of 
plasticity, in which the plastic strain increment is 
determined by the stress and the stress increment. Rate- 
dependent plasticity models proposed in the postwar period 
included those by Freudenthal [34], Malvern [35,36]* who 
constructed a one-dimensional equation relating the stress 
increment to the strain rate, Lubliner [37], who modified 
Malvern's equation to include a limiting stress-strain 
curve, and Perzyna [ 38 - 41], who generalized Malvern's 
equations to three dimensions and finite strains. 

Many classical plasticity models have been extended to 
include rate dependence by adding an uncoupled rate- 
dependent strain component to the rate-independent strain 
component. No observable difference exists between the 
physical mechanisms causing these phenomena, however, and 
it has been argued that the two terras should be combined 
[42]. Naghdi and Murch [43] attempted to include rate 
dependence by introducing interdependent, although not 
unified, plastic and viscoelastic strain components. 

Classical plasticity theories are usually described by 
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i jmn 
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mn 
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(32) 


T 

where D. . is the linear elastic modulus tensor, e 

ijran mn 

the thermal strain, and e 1 , the inelastic strain, is 

ran 

given by 


i s 


• I 
e mn 


A 


oF/o o 


ran ' 


(33) 


where A is a scalar valued function of state and F is a 
scalar valued function for Inelastic behavior often taken 
to be the yield function. If F is described by the von 
Mises yield criterion [2*0, given by 


FU lj * “lj ) “ 2 ( 0 1 j ' “lj )( ”lj • “lj 1 ' k2> <31|) 

where is a tensor describing the yield surface center 

in stress space and k is a constant representing the yield 
surface size, then equations (33) can be written as 



A 




(35) 


resulting in a kinematic hardening model with constant 
yield surface size. Furthermore, if the yield surface 
translation is derived from the Ziegler modification [32] 
of the Prager work-hardening rule [31 ]» it may be described 
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by 


a 


■tj 




(36) 


where p is a scalar valued function of state. 

By use of equations (32), equations (35) and (36) can 
be shown to be consistent with growth laws (15). It can 
then be seen that classical plasticity theories are 
consistent with the thermodynamic framework developed in 
the equation of state approach. 

Recently proposed theories include those by Snyder and 
Bathe [44], Yamada and Sakurai [45], Allen and Haisler 
[46,47], which attempts to model transient temperature 
effects, Zienkiewicz and Cormeau [48], and Robinson [49]. 
Bodner and his associates proposed a model [50] which uses 
a flow law similar to equations (31), although the internal 
state variable growth laws are based on microstr uc tural 
considerations. Robinson's model and Bodner's model are 
still under active development. 


Nonlinear Viscoelasticity Models 
Nonlinear viscoelasticity models are usually based on 
thermodynamics and/or mechanical analogs composed of 
springs and dashpot combinations such as the three 
parameter standard solid, shown in Figure 4. In 1954, Biot 
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[6,7] proposed a theory of linear thermoviscoelasticity 
based on the principles of irreversible thermodynamics. 
Schapery [12] later modified this to include nonlinear 
behavior and Introduced a reduced time scale to solve the 
resulting highly nonlinear equations. Valanis [51-53] 
proposed a viscoplasticity theory, termed the endochronlc 
theory, which includes a history-dependent reduced time 
parameter. This was later shown to be a special case of 
Schapery's model [54]. Some other models based on 
thermodynamics are those of Coleman and Noll [55,14], Green 
and Rivlin [56], and Green and Naghdi [57]. These tend to 
be quite complex in nature and difficult to solve. 

Many of the recent models are based on variations of 
the three parameter standard solid. The standard solid 
shown in Figure 4 has been demonstrated to be an 
appropriate representat ion for therinov i scoplast ic metals if 
the springs and dashpot are nonl i near i zed [58,59]. The 
governing differential equation for this solid is 

0 + (n M /E M } ° = E ® £ + (n M /E M } (E M + E oo )£ ' (37) 


This may be rewritten in the equivalent form 


= o / ( E + E 
M 


• } + e m /u m (e m 


♦ E_ ) ] ( o 


E e) . (38) 

CO 


e 



In accordance with the instantaneous linear behavior of 
metals, it is assumed that 


E w *■ E = E = Elastic modulus =» constant, (39) 

M ® 


and equation (38) becomes 


(40) 


where e 1 is the inelastic strain, given by 


e = (E M /n M E) ( o - E e). (41) 

MM ® 


The nonlinearity of the springs and dashpot can be 
introduced by considering a multiaxial extension of 
equat ion (37): 


pq 


^pqmn°mn 


G e + 
pqmn mn 


M e _ , 
pqmn mn 


(42) 


In order to model metals, K , G , and M are 

pqmn pqmn pqmn 

nonlinear material functions determined from experimental 
data. In addition, constraint (39) requires that 


K ijmn M mnkl D ijkl’ 


(43) 
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where ^ la the linear elastic modulus tensor. 

Equations (42) and (43) lead to 


I T 

o . . = D , . (e - e - e 
ij ijran ran ran mn 


where the inelastic strain e 


ij 


e ij " M ijkl (o kl ' G klmn e mn ) 


) , 

is given by 


(44) 


( 45) 


Substitution of equations (44) into (45) will result in 
equations consistent with growth laws (15). It can then be 
seen that nonlinear viscoelasticity models are consistent 
with the thermodynamic framework developed earlier in the 
equation of state approach. 

Some recent models which are based on nonlinear 
modifications of the three parameter standard solid are 
those of Krempl, et al . [59-66], and Walker [58,67,68]. 


M icrostr uctural Models 

Microstr uctural models have generated increased 
interest in recent years because their formulation is based 
on microphysical considerations rather than a statistically 
averaged macroscopic view. The primary micromechanisms 
which cause inelastic deformation are dislocation glide and 
climb, deformation twinning, diffusion, grain boundary 
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sliding, and microvoid growth and coalescence. These 
mechanisms are difficult to model individually and are 
usually reduced to three more general mechanisms: 
dislocation arrangement, dislocation density, and damage. 

In the concept of microphenoraenology it is assumed that the 
three mechanisms can be character ized by internal state 
variables. One internal state variable represents the back 
stress for modelling the dislocation arrangement. Another 
internal state variable, termed the drag stress, models 
dislocation density. In addition, it can be argued that 
the inelastic strain represents the third internal state 
variable [15,16,19]. Current research is in progress to 
account for damage in materials using a fourth internal 
state variable [69-71], although it is difficult to 
distinguish between damage and the drag stress since both 
are primarily stiffness-reducing parameters. 

Early research in this field was performed by Coble 
[72], Nabarro [73], and Herring [74] in their work on 
d i f f u s i on- con tr ol 1 ed creep. Other important contributions 
were made by Sherby, et al . [75-77], Garofalo [78], Argon 
[42], Weertman [79], Alden [80], Kocks [81], and Hart [82]. 

Mukherjee, et al . [83], studied the effectiveness of 
using a power law of the following form to model 


dislocation climb: 


where D is the sel f -d i f f us i v 1 t y , E is Young’s modulus, and 
A and n are temperature-dependent material parameters. The 
authors concluded, however, that a model in terms of a 
single dislocation mechanism may not completely 
characterize high- temperature creep. 

Gibbs [84] proposed a two- i nt ernal- s ta te- var iable 
theory in which the inelastic strain rate is characterized 
by an exponential function of temperature and an effective 
stress: 

e 1 = A* exp{-[Q c - ( o- o a ) 1 b 2 ] / kT } , (47) 

* 

where A is a structure factor, Q is the dislocation core 

c 

diffusion energy, k is the activation energy, T is the 

absolute temperature, b is the Burger's vector, and o and 

1 are internal state variables representing, respectively, 

back stress and drag stress. 

It is generally recognized [85] that at intermediate to 

high stresses and at temperatures above 0.5T , where T is 

mm 

the absolute melting temperature, stress is related to the 
inelastic strain rate by a power law of the form 
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where A and n are constants. At very high stress levels, 
the relationship is 

e 1 = B exp ( mo ) , (49) 

where B and m are constants. Power law and exponential 
relationships represent limiting cases for the more general 
empirical relationship 

. i ' ' 

\ 

e 1 - C [sinh(k 0 )] P , (50) 

where C, p, and k are constants. Equation (50) 
approximates a power law when ko<0.8, but it approximates 
an exponential curve when ko>1.2. See Figure 5 for a 
comparison of the forms of equations (48), (49), and (50). 

Microstructural models are usually similar in nature, 
differing only in the proposed relationship between stress 
and inelastic strain (that is, equations (48), (49), or 
(50)) and in the proposed internal state variable growth 
laws. The growth laws for the internal state variables 
representing the back stress and the drag stress usually 
follow the framework established by Bailey [86] and Orowan 
[87] in which a hardening term, proceeding with accumulated 
deformation, competes simultaneously with a softening or 
recovery term, proceeding with time. For example, the 
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isothermal model proposed by Krieg, et al . [88], is given 
in its uniaxial form by 

o = E (e - e 1 ) , (51 ) 

where the inelastic strain is given by 

e 1 *■ c 1 [ | o - a.||/a 2 ] sgn(o - a 1 ) , (52) 

where the growth laws for the back stress a 1 and the drag 
stress a 2 are given as the differences of hardening and 
softening rates: 

* * I 2 2 - 

a 1 = c^e - c iJ a 1 [exp(c 5 ap - 1] sgn(a 1 ) (53) 

and 

a 2 = c 6 l eI | ” °7 ^ a 2 ~ a 2 ^ ’ ( 54 ) 

where c 1 , c 2> c^, c^, c^, , c? , n, and a 2 are 

temperature-dependent material constants, and sgn() is the 
signum function. 

Substituting equation (51) into (52) will result in 
equations consistent with growth laws (15). In addition, 
since equations (53) and (54) are consistent with growth 
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laws (15), It can be seen that rnicrostructural models agree 
with the thermodynamic framework developed earlier in the 
equation of state approach. Cescotto and Leckie [89] 
proposed a model which has the interesting feature that no 
particular forms are assumed for the inelastic strain rate 
function and the hardening and recovery functions; only 
weak hypotheses are required and the functions are defined 
experimentally. 

Other recently proposed models include those of Miller, 
e t al . [90], in which a hyperbolic sine function is used 
for the stress- inelast ic strain rate equation, Walker [58], 
in which a power law is used for the inelastic strain rate 
equation, and Bodner [50], in which an exponential function 
is used. An interesting aspect of these last two models is 
that, although the respective growth laws for the back 
stress and the drag stress are formulated from 
rnicrostructural considerations, Walker's inelastic strain 
rate equation is based on a nonlinear modification of a 
three parameter viscoelastic solid, while Bodner's 
inelastic strain rate equation is based on a flow law 
similar to the Prandtl-Reuss equation of classical 
plasticity (see the literature reviews on nonlinear 
viscoelasticity models and classical plasticity models, 
respectively) . 


Since theories based on rnicrostructural cons i der a t i ons 
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have been given a great deal of study in recent years, the 
theories chosen for comparison are taken from this field 
and are those of Krieg, et al . [88], Miller, et al . [90], 
and ,Bodner, et al . [50]. The motivation for the choice of 
these three models is that they are still under 
development; they have been cast in a common thermodynamic 
framework; the authors have included in their theories 
detailed accounts of methods of determination of the 
material constants; these models have received considerable 
attention [20,68,91-93]; and research is in progress to 
extend these isothermal models to include transient 
temperature response [68,91.94]. 
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KRIEG, SWEARENGEN, AND ROHDE’S MODEL 

Chronological Development 
In 1977, Krieg [95] cast several current unified 
constitutive models into the common framework of a skeletal 
model in which general similarities, such as the use of 
competing hardening and recovery rates for internal state 
variable growth laws, are noted. These models include 
those of Bodner, et al . , Robinson, et al . , Hart, Lagneborg, 
Miller, Pasley and Wells, Ponter and Leckie, and an early 
development of Krieg's model, which is described as a 
composite of the aforementioned models. All hardening and 
softening functions are taken to be constants in the 
skeletal model in order to make observations about 
difficulties encountered ii numerically integrating 
mathematically stiff constitutive equations. Finally, 
integration techniques are proposed to alleviate these 
diff iculties. 

In a development of this model, Krieg, Swearengen, and 
Rohde [88] proposed a power-law kinetic relation in 
multiaxial form to model pure aluminum at room temperature. 
Th.s model has three internal state variables: the 
inelastic strain, the back stress and the drag stress. The 
growth laws of the back stress and the drag stress are 
formulated in the usual hardening rate/recovery rate 


format. The hardening functions have not been defined and 


are assumed to be constant; recovery is considered to be a 
thermally-activated process and dynamic recovery is not 
included. Furthermore, the response of pure aluminum at 
room temperature is very nearly kinematic and, as a 
consequence, the drag stress (the isotropic hardening 
variable) is taken to be constant. Krieg, et al . , present 
a test method in which all material constants, except 
Young's modulus, are calculated by determining values of 
the back stress from a set of stress-drop tests. In this 
test, a portion of the load on a specimen during secondary 
creep is removed and the resulting strain transient is 
recorded. The authors have obtained reasonable results for 
pure aluminum at room temperature for various load 
histories, although the predicted hysteresis loops are 
"over-square" and the model's use is limited in high strain 
rate ranges. 

In 1982, Jones, et al . [96], noted that, while most 
unified creep-plasticity models require an extensive number 
of tests for characterization of material parameters, the 
model proposed by Krieg, e t a 1 . , requires only a set of 
stress-drop tests. It is pointed out that, although simple 
in concept, the stress-drop test is difficult to perform in 
practice. The critical measurement in this test, a zero 
strain rate after unloading or merely the absence of a 



35 


resolvable strain rate, results in a strong dependence of 
the parameters on experimental resolution. The authors 
present the design of a test apparatus to minimize these 
effects. In addition, there appear to be differing 
opinions in the literature on the response of the creep 
rate immediately after unloading. Jones, et al . , conclude 
that stress-drop testing should be performed with great 
car e . 

In subsequent research in 1983, Jones and Rohde [97] 
refer to a new technique to reduce stress-drop data which 
alleviates the difficulties mentioned above. Rather than 
attempting to measure a zero strain rate, this method [98] 
uses a comparison of the instantaneous strain change and 
the maximum strain change after unloading. This method 
appears to simplify the calculations considerably and 
merits further study. 

General Theory 

In their unified creep-plasticity model, Krieg, 
Swearengen, and Rohde [88] state that inelastic deformation 
can be caused by: 1) thermally-activated, stress-assisted 
penetration of short-range obstacles to dislocation motion, 
2) generation and immobilization of dislocations, leading 
to increases in dislocation density, stored elastic energy, 
and the flow stress, and 3) stress-assisted thermal 
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rearrangement of dislocations to reduce long-range internal 
stresses. Krieg, et al . , note that strain hardening at low 
to intermediate homologous temperatures (below 0 . 51 ^, where 
T m is the absolute melting temperature) is both isotropic 
and kinematic. Isotropic hardening, manifested by an 
increase in the height of the reversed strain hysteresis 
loop, has its microstructural origins in dislocation 
obstacle interaction, while kinematic hardening, manifested 
by the Bauschinger effect, has its origins in dislocation 
pile-ups. Mathematical descriptions of this behavior 
require two internal state variables. 

Deformation kinetics resulting from thermally-activated 
processes are usually described by Arrhenius or hyperbolic 
sine functions; however, over a wide range of stresses, the 
analytically simpler power law provides a good 
approximation . This motivated Krieg, e't al . , to propose the 
following flow rule in multlaxial form in deviatoric 
stress-strain space: 

S 1 " e 0 < I 5 I /R > m I » ( 55 ) 

. i , , 

where e is the inelastic strain rate vector, || is the 
Euclidean vector norm, and m are temperature-dependent 
material constants, and R is the drag stress, an internal 
state variable which accounts for isotropic hardening. The 
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effective stress £ is given by 


5 = 3 


a , 


( 56 ) 


where a is the back stress, an internal state variable which 
accounts for kinematic hardening; and s is the applied 
stress. The increase in flow stress associated with an 
increase in hardening is taken to be proportional to inelastic 
strain. Recovery is associated with dislocation escape from 
glide planes by climb or cross-slip. In the absence of 
external loads, this escape reduces the internal stresses 
through annealing; in the presence of external loads, this 
escape permits additional inelastic deformation through creep 
or dynamic recovery. The evolution of a and R is given by 
the difference between hardening and recovery rates as 


a = A e* - r ( <*/ | 3 | ) 


(57 ) 


and 


R - A R |e - r R , 


( 58 ) 


where A and A n are the hardening functions, and r and r D 
a it an 

are the recovery functions. 


Krieg, et al . , assume the hardening functions A q and Aj 
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to be constant. The recovery rates are assumed to be 
proportional to dislocation velocity of escape and 
inversely proportional to the escape distance [99], where 
the velocity and distance depend on the magnitude of the 
back and drag stresses. Since annealing rate depends on 
hardening mechanisms, one would expect the kinetics of 
recovery for a and for R to be different. A unique 
feature of this model is that it incorporates different 
physically-based recovery kinetics for each internal 
variable . 

Dislocation processes associated with the drag stress 
are formations of dislocation tangles, or networks. 
Temperature-dependent changes in the networks can be 
described by the climb recovery model of Friedel [99]: 

r R = -K 1 ( R n /kT ) exp(-U/kT) , (59) 

where and n are material constants, k is the gas 
constant, T is the absolute temperature, and U is the 
activation energy associated with the particular 
micromechanism of softening. 

Dislocation processes associated with the back stress 
are pile-ups or cell-wall bowing. Although both cross-slip 
and climb are recovery mechanisms, only screw dislocations 
can cross-slip and, consequently, the model's authors feel 
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that recovery by climb is a more complete description. The 
kinetics of this process are again given by Friedel [99]: 

r^ = -K 2 a 2 [exp(K 3 a 2 /kT) - 1], (60) 

where and are material constants. 

By U3e of equations (60), in multiaxial form, and (59), 
equations (57) and (58) become 

a = A a g I - | a j K 2 [exp(K 3 | a 2 | /kT) - 1] (61) 

and 

R = A R | I 1 | - K 1 [(R - R Q ) n /T] exp(-U/kT), (62) 

where Rq represents an isotropic annealed state. 

For conditions of uniaxial stress and constant 
temperature, equations (55), (61), and (62) reduce to 

e 1 = c 1 j ^ / R j ° 2 sgn(c), (63) 

a = - c^a [exp(c^a ) - 1] sgn(a), (64) 


and 
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R = c 6 j c 1 j - c ? (R - R 0 ) n , (65) 

where , c 2> c^, c^, e^, , c^ , n, and R Q are 

temperature-dependent material constants, sgn () is the 
signum function, and the effective stress q is 

5 •• o - a . ( 66 ) 

Krieg, et al . , simplify the model further by assuming a 
kinematic work-hardening material so that the isotropic 
hardening variable R will remain constant. Equation (65) 
can then be omitted and the uniaxial isothermal form of the 
model becomes 

o = E ( e - e 1 ) , (67) 

where e 1 is given by 

e 1 = Cl | c| c 2 sgn ( 5 ) , ( 68 ) 

and 

• • X 2 2 

a = c^e " c 4 a [exp(c,.a ) - 1] sgn(a), (69) 

where c, , c_, c_, c „ , and c c are-temperature dependent 

1 ^ 3 M o 



in 


material constants, and the effective stress q is given by 
equation (66). 


Evaluation Of Material Constants 
Krieg, et al . , present an experimental test and data 
reduction procedure which allows all constants, except the 
elastic modulus E, to be determined from a set of 
stress-drop tests, sometimes called strain transient dip 
tests . 

Before describing the procedure to evaluate the 
material constants, it would be instructive to review the 


somewhat controversial stress-drop test. High-temperature 
creep can be characterized in terms of an internal stress 
and an effective stress [1003, where the internal (or back) 
stress is the driving force for recovery, and the effective 
stress (the difference between the applied and back stress) 


is the driving force for dislocation glide. 


Consequently, 


it is of interest to be able to determine the back stress. 


Since the back stress is a result of the dynamic 
balance between strain hardening and recovery, the back 
stress begins to change when the applied stress is changed. 
Therefore, the transient response examined should be of 
short enough duration to ensure a structure which is still 
reasonably representative of the steady-state structure; 
i.e., a structure in which the value of the back stress has 
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not changed appreciably from its steady-state value. A 
technique described by Ahlquist and Nix [100,101] involves 
the dynamic measurement of back stresses by rapidly 
reducing the applied stress from the steady-state value 
and measuring the strain rate immediately after the 
reduction. For a small reduction in stress Ao, the strain 
rate immediately after the reduction is positive as shown 
in Figure 6a. For a large reduction in stress Ao, the 
strain rate is negative as shown in Figure 6b. For an 
intermediate reduction in stress Ao, the strain rate is 
zero and the new stress level o 2 describes the mean 
internal (or back) stress o^. After a short time, recovery 
events reduce the back stress, yielding a positive strain 
rate. As can be seen from Figure 7, it may be necessary to 
load and unload a sample several times during a test before 
finding the stress reduction that will give the back 
stress. 

This procedure for determining the back stress applies 
to thermally-activated creep. However, there are different 

i 

hypotheses for creep deformation where the concept of a 
back stress may become questionable. Consequently, these 
hypotheses give differing predictions of the transient 
response after a stress drop. Poirier [102] summarizes 
these into four categories, as shown in Figure 8, where the 


responses due to stress drops of different magnitudes have 








been superposed. 

The athermal recovery-controlled creep case is 
characterized by: 1) the absence of an effective stress, 

2) creep proceeding only if the back stress is reduced to 
the applied stress level by diffusion-controlled recovery, 
and 3) stress reductions, regardless of magnitude, always 
followed by a period of zero creep rate (see Figure 8a). 

The thermally-activated creep, jerky glide, case is 
characterized by localized obstacles to dislocation 
movement and, consequently, no negative creep rate after 
the stress reduction. However, there may be a period of 
zero creep or positive creep, depending on the magnitude of 
the stress reduction (see Figure 8b). 

The thermally-activated creep, viscous glide, case is 
characterized by "smeared" rather than localized obstacles. 
Consequently, there may be positive, zero, or negative 
creep after the stress reduction (see Figure 8c). 

Finally, the thermally-activated creep case with a 
combination of jerky and viscous glide is characterized by 
a combination of localized and smeared obstacles. The 
response of this case is similar to that of the previous 
two and depends on the magnitude of the stress reduction as 
well as the applied stress level (see Figure 8d). 

Poirier concludes that, while the concept of a back 
stress cannot be rejected on experimental grounds, 
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Fig. 8 (continued) 
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measurements of the back stress resulting from stress-drop 
tests should be performed with great care and are valid 
only for thermally-activated flow. 

Similarly, Jones, et al . [96], conclude that due to 
ambiguity involved in attempting to determine a zero creep 
rate or merely the absence of a resolvable creep rate after 
unloading, the stress-drop test should be Interpreted with 
great care on the basis of a specific deformation 
hypothesis. With this in mind, one can now proceed to 

evaluate the material constants c, , c„, c_ , c,., and c r . 

1 c. 5 *4 b 

c 1 and c 2 

Equation (68), which governs inelastic strain, can be 
rewritten as 

In ( e * ) = ln(c 1 ) + c 2 ln|cj, (70) 

where |c| is the difference between the applied stress o 

. i 

and the experimentally determined back stress a, and e is 
the creep rate immediately preceding the particular stress 
reduction associated with that value of a. The left hand 
term of (70) must be linear in ln|c| if it is an adequate 
representation of experimental data. A plot of data in the 
form of (70) is shown in Figure 9. Graphical means or a 
least-squares analysis of the data will give the values of 
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c ^ and • 

and (c^/Cjj) 

. T 

During steady-state creep, the creep rate e is 
constant and, as a consequence of equation (68), £ must 
also be constant. Since the applied stress o is constant 
in a creep test, this Implies that the back stress a is 
also constant. With a = 0, equation (69) reduces to: 

(c 3 /c J| ) = La 2 / ( c 1 | c | c 2 ) ] [exp ( c^a 2 ) - 1], (71) 

where the fitted form c ^ ^ | 2 has been used for the term 

e 1 to reduce experimental scatter. The experimental pairs 

(^, 0 ^) and ( ^ 2 ' a 2 ) ' where the subscripts 1 and 2 refer to 

the highest and lowest values of a, may be used in equation 

(71) to give two nonlinear equations in two unknowns. 

These are then solved by iterative means to give the values 

of c,- and the ratio (c_/c„). 

5 3 4 

and c^ 

It now remains to determine either c^ or c^, since only 
their ratio is known at this point. In order to accomplish 
this, information from a point on the primary creep curve 
(where a is not zero) must be used since the available 
steady-state information has been exhausted. Substituting 
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equation (68) into (69), separating variables, and 
integrating leads to 

c 3 fc “ ; 0 l {c i I °i " a l° 2 “ 

(c 1< /c 3 )a 2 [exp(c 5 a 2 ) - I]}' 1 da, (72) 

wher’e is the creep stress, a is the variable of 
integration, and t are the back stress and time, 

. i 

respectively, at which the primary creep rate is still 
twice the eventual steady-state value. The value for at 
this data point can be determined by inverting equation 
( 68 ) : 

- o, - (e 1 I /c 1 ) W ° 2 . (73) 

Since all the quantities inside the integral now have known 
values, equation (72) can be integrated numerically to give 
a value for c^t. Since t is known, c^ and, consequently, 
c ^ can be calculated easily. Although not explained by 
KrieB, et al . , a primary-to-steady-state creep rate ratio 
of two appears to represent a transition point between the 
region of rapidly decreasing primary creep rate and the 
region in which the creep rate slowly approaches its 
steady-state value, as shown in Figure 10. 



Fig. 10 A qualitative graph of the pr i mary- to- secondary 
creep ratio used in determining c_ and c^ 
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In summary, the Krieg, Swearengen, and Rohde theory 

requires a constant strain rate test (for the elastic 

modulus E), several stress-drop tests in the steady-state 

creep region (constants c, , c„, c c , and the ratio c./c.,) 

\ d o 3 4 

and a creep curve complete with primary and secondary 
regions (constants c^ and c^). This concludes the 
determination of the material constants for the Krieg, 
Swearengen, and Rohde model. 



BODNER AND PARTOM'S MODEL 


Chronological Development 

In 1975, Bodner and Partom [50] expanded an earlier 
model, capable of modelling e 1 as t i c- per f e c tl y plastic 
behavior only [103], to include strai n-hardening by 
introducing plastic work as the measure of the hardened 
state. The formulation of the inelastic strain rate 
equation is based on a generalization of the flow rule of 
classical plasticity and is motivated partly by dislocation 
dynamics. The inelastic strain rate is proposed as an 
exponential function of stress and an internal state 
variable representing hardness, where the hardness 
parameter is an exponential function of plastic work. The 
proposed formulation corresponds to isotropic hardening and 
therefore would not characterize the Bauschinger effect. 
Bodner and Partom found good agreement of their theory with 
experimental results for commercially pure aluminum. 

In noting that most constitutive theories consider 
isotropic and directional hardening effects as completely 
separable and thereby controlled by different internal 
state variables, Bodner and his colleagues [104] extended 
their isotropic theory in 1979 to model uniaxial cyclic 
loading by proposing an alternative approach. This 
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approach uses separate values of the hardening parameter in 
order to account for the directional character of hardening 
during cyclic loading; one value corresponds to tensile 
loading, the other corresponds to compressive loading. 
Bodner, et al . , found good agreement with uniaxial 
experimental data for OFHC copper and commercially pure 
titanium and aluminum at room temperature in modelling 
phenomena such as cyclic strain hardening and softening. 

The Bodner- Par tom theory subsequently was generalized 
to a framework for the multidimensional stress case by 
using a general anisotropic form of the flow law [105]. 

In another development of the theory [106], an 
additional term was introduced to the constitutive 
equations in order to account for the thermal recovery of 
hardening, a term essential for modelling secondary creep. 
This development was used to model the superalloy Rene 95 
at high temperatures. 

A parameter to account for damage was introduced in 
1980 [71]. Bodner examined the isotropic case first and 
then generalized the damage parameter to the anisotropic 
case in a manner similar to his earlier treatment of 
anisotropic hardening [105]. 

A preliminary method to determine the material 
constants of the uniaxial constitutive equations directly 
from experimental data rather than from tr ial-and-error 
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c ur ve - f i t t i ng was presented in 1981 in a study of Rene 95 
at high temperatures [69]. This was further developed into 
a detailed systematic method and applied to IN100 at 
elevated temperatures [107] for the case of isotropic 
hardening with thermal recovery and no damage. 

In an earlier paper [105], anisotropy was introduced in 
a manner not automatically leading to plastic 
incompressibility. Bodner and Stouffer revised the theory 
in 1983 [108] to enforce plastic incompressibility in order 
to make the theory consistent with stability and 
thermodynamic principles. Results from the revised theory 
are essentially unchanged from the original work except 
that plastic volume changes become zero. 

An incrementally isotropic form of the flow law has 
been proposed [109] in order to alleviate some of the 
computational difficulties encountered in using the full 
anisotropic form. In this formulation, the scalar 
hardening variable is taken to be the sum of the isotropic 
component and a scalar effective value of the directional 
component. Evolutionary equations are given for both 
components, including thermal recovery terms, as well as 
for isotropic damage. It is suggested that anisotropic 
damage can be treated in a manner similar to anisotropic 
hardening. 


Recent work on the Bodner-Partom theory includes 
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research into modelling transient temperatures [91, 94], a 
complete detailed review of the present constitutive theory 
[110], and development of a systematic method for the 
determination of material constants for the model in its 
complete anisotropic form with damage [ill]. 


General Theory 

The unified creep-plasticity model of Bodner and his 
colleagues [107] is motivated in part by work done in 
dislocation dynamics in which the total deformation can be 
separated into elastic and plastic components which are 
generally non-zero at all stages of loading. A yield 
criterion or loading and unloading conditions are therefore 
not required. 

The total strain rate is given by 


e 


i J 



(7 4 ) 


• E 

where the elastic strain rate e_ is given by the time 
derivative of Hooke’s Law and the inelastic strain rate 
is given by the flow rule of classical plasticity: 


• I -1 . 

= e.. = A o!., 

A J A J A J 


(75) 


where e. and o', are the deviatoric components of the 
a J a J 
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inelastic strain rate and applied stress tensors, 
respectively, and A is a scalar material function. 
Squaring both sides of (75) gives 


A 2 = D 2 /J 2 , (76) 

p 

where D is the second invariant of the inelastic strain 
rate deviator 


D 


P 

2 


1*1 • I 
2 e ij e ij 


(77) 


and J 2 is the second invariant of the stress deviator 


f 2°ij°ij* 


(78) 


A fundamental assumption of the Bodner-Partom theory is 
that all inelastic deformations are governed by the kinetic 
equat ion 


°2 ' f < J 2- T 'V’ 


(79) 


where T is the absolute temperature and Z ^ are a set of 
internal state variables. Extensive work in the field of 
dislocation dynamics [11 2-1 I 1 *] has shown that dislocation 
velocity and, therefore, inelastic strain rate can be 
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represented as a power function or an exponential function 
of stress. 

Bodner and Partom note that in order to represent 

p 

material behavior should vary inversely with the 
measure of strain hardening, have a limiting value for 
large and be almost zero for low stresses. Whereas 

several mathematical functions may fulfill these 
requirements, the function chosen by Bodner and Partom is 

D 2 = D 2 exp[-(Z 2 /3J 2 ) n ] , ( 80 ) 

where D Q is the limiting strain rate in shear, n is a 
temperature dependent material constant, and Z is 
interpreted as an internal state variable representing 
hardness. It should be noted that, in subsequent 
developments of the theory, an internal state variable 
representing damage was introduced in the following manner 
[71,109,110]: 

D P 2 = D 2 exp { - [ Z 2 ( 1 - u)) 2 /3J 2 ] n l, (81) 

where w represents damage in the material. 

Equations (75 ), (76 ), and (80) can then'be combined to 
give an expression for the inelastic strain rate: 
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i\. - {D 2 exp[-(Z 2 /3J 2 ) n ]/J 2 } 1 /2 o [y (82) 

The work-hardened state, that is, the resistance to 

plastic flow, is assumed to be represented by a single 

P P 

variable Z which depends on the plastic work W , where dW 
is ode*. The plastic work rate was chosen by Bodner and 
Partom as the measure of hardening over the more common 
choice of the inelastic strain rate because Z is directly 
related to the stored energy of cold work and this choice 
leads to relative simplicity of the resulting equations 
[ 110 ]. 

It is postulated that the evolution of Z is governed by 
the current values of stress, hardness, and the absolute 
temperature : 

Z = F( J 2 ,Z,T) , (83) 

where the particular form chosen by Bodner and his 
colleagues is 

Z = m (Z 1 - Z) W P - A Z [(Z - Z 2 )/Z 1 ] P , (84) 

where m, , Z ^ , A, and r are temperature-dependent 
material constants and the constant Z^ is designated as the 
initial value of Z. The constants Z^ and Z ? correspond to 
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the maximum value of Z and the minimum recoverable value of 
Z, respectively. 

The first term in equation (84) i3 the hardening term, 
where the negative part can be interpreted as dynamic 
recovery, and the second term is the thermal recovery term, 
designed to be negligible during rapid loading histories. 

This representation, equation (84), corresponds to 
isotropic hardening and, as a result does not predict the 
Bauschinger effect. The Bodner-Partom theory was later 
modified [109] to include directional hardening by taking 
the scalar hardening variable Z in the inelastic strain 
rate equation (82) as the sum of the scalar, isotropic 
hardening variable Z 1 and a scalar effective directional 

A 

hardening variable Z : 

Z = Z 1 ♦ Z A , (85) 

where Z 1 is now given by (84) and the evolution equation 

A 

for Z has the general form of (84) but is tensorial in 
character. 

For the purposes of this research, the isotropic 
hardening model with no damage [107], equations (82) and 
(84), is considered. 

For uniaxial stress and constant temperature, equation 


(82) reduces to 
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c 1 = D q exp[-^(o/Z) 2n j sgn(o) (86) 

and equation (84) remains unchanged due to its scalar 
nature : 

Z = m (Z - Z) W P - A Z [( Z - Z 2 )/Z 1 ] r . (84) 

The material constants to be determined are , n, m, Z 1 , 
A, Z , r, and Z Q , the initial value of Z. 


Evaluation of Material Parameters 


The quantity D Q is the maximum value of strain rate in 
shear and has its physical basis in the upper bound of the 
dislocation velocity. The value of D Q can be set as 

1 0 8 sec -1 for e I >10 3 , 10 6 sec" 1 for 10 1 <e I <10 3 , and 

4 - 1 • I 1 

1 0 sec for e <1 0 . 


Z and n 

The first step in the evaluation of material constants 
is to determine the strain rate sensitivity parameter n and 
the maximum value of the hardness, Z , from tensile data. 
For material behavior as shown by curve (a) in Figure 11 
where the stress saturates to a maximum value, in the 


region where both the stress and strain rate are constant Z 
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must also be constant in order to satisfy equation (86). 

For short duration tests with no recovery, the material 
must be in a fully work-hardened state to obtain the 
maximum value of stress; that is, Z must have its maximum 
value Z . Rewriting equation (86) for this steady flow 
condition gives 

ln[-ln(^|e I /D 0 )] = 

-2nln(o) + [2nln(Z ) ♦ ln(|)]. (87) 

The left hand term of (87) must be linear in In (o) if it 
is an adequate representation of experimental data. A plot 
of data in the form of equation (87) is shown in Figure 12. 
Graphical means or a least-squares analysis of the data 
will give the values for n and Z^ . 

However, if the tensile curves are not quite saturated, 
as shown by curve (b) in Figure 11, simply using the 
largest value of stress attained may give incorrect values 
for Z^ and n. In this case, equation (86) is solved for o 
to give 

o = K 1 Z, (88) 


w lie re 
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K 1 - [21n( 7 |D 0 /e I )]“ 1 /2n . (89) 

In the absence of hardening recovery, the hardening 
evolution equation (8*0 becomes the first order linear 
differential equation 

dZ - m (Z 1 - Z) dW P , - (90) 

which can be combined with the differential of equation 
( 88 ) to give 

do = K 1 m ( Z - Z ) dW P . (91 ) 

By defining a new parameter Y as 

Y = do/dW P (92) 

and substituting equation (88) for Z, equation (91) can be 
rewritten as 

Y = K m Z - m o. (93) 

In order to be an adequate representation of experimental 
data, y must be linear in o, as shown in Figure 13* The 
saturation stress is given by the stress level where Y=0. 
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With this value of the saturation stress, one can now 
return to the method outlined above to determine n and . 

Z and m 

The value of the hardening variable Z can be calculated 
for corresponding values of stress and inelastic strain 
rate by inverting equation (86): 

Z = o C21n( 7 |D () /e 1 )] 1 /2n . (94) 

Equation (90), the hardening evolution equation with no 
recovery, can be integrated to give 

ln(Z 1 - Z) - ln(Z 1 - Z Q ) - mW P , (95) 

where Z Q is the initial value of Z. Several values of Z 
can be calculated for each test by using equation (94) at 
different stress levels and their corresponding values of 

. i 

e taken from the region in the a~c graph where recovery 

p 

is minimal. In this manner a plot of ln(Z - Z) against W 
can be made for several constant strain rate tensile tests. 
Since this plot must be linear in order to represent 
material behavior adequately, as shown in Figure 14, the 
constants m and Z ^ can be calculated easily either by 
graphical means or by a least-squares analysis. 
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Alternatively, if a Y - o graph was used in the earlier 
determination of n and , then the average of the slopes 
of the linear curve-fits of the Y-o data will be -m. Z Q is 
then given by calculating values of Z in the immediate 
post-yield (i.e., post-elastic limit) region of a 
stress-strain curve. 

Z 2 

During secondary creep, since stress and inelastic 
strain rate are constant, the value of Z must also be 
constant in order to satisfy equation (86). These 
stationary values of Z can be determined from equation (9*0 
for several creep tests. The lowest calculated stationary 
value of Z can be used for the value of Z 2> the minimum 
value of Z, although Bodner states that Z 2 = Z Q for many 
applications since the initial condition (Z Q ) is the fully 
recovered state (Z^) [110]. 

A and r 

Since Z = 0 during secondary creep, the hardening 
evolution equation becomes 

. m (Z - z) W p = A Z 1 [(Z - Z 2 )/Z 1 ] r . (96) 

. i 

The values of o, e , and Z used in the determination of Z^ 
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• P 

above may be used in a plot of ln[m ( Z - Z) W ] against 
ln[(Z - Z 2 )/ Zl ]. The constants A and r can be calculated 
from these data by assuming a linear representation, 
although a bilinear representation may also be possible in 
some case s [ 1 07 ] • 

In summary, the Bodner- Par tom theory requires several 
constant strain rate tensile tests (constants E, n, Z^ , m, 
and Z Q ) and several creep tests (constants Z^ , A, and r), 
where the actual number of tests is dependent on the amount 
of data deemed necessary for accurate evaluation of the 
constants. This concludes the determination of material 


constants for the Bodner- Par tom theory. 



MILLER'S MODEL 


Chronological Development 

Miller's model [90], introduced in 1976, is based on a 
combination of creep phenomenology and microstructural 
considerations. In stressing the need to model 
steady-state creep rates accurately, Miller proposes to use 
the Garofalo hyperbolic sine relation for steady-state 
creep as a basis for his constitutive equations. 
Representing isotropic and kinematic hardening, 
respectively, a drag stress variable and a back stress 
variable are introduced to describe non-steady-state 
inelastic responses. The growth laws for the back stress 
and the drag stress follow the standard Bailey-Orowan 
format of competing hardening and recovery rates. Miller's 
growth laws feature a constant kinematic work-hardening 
coefficient but a variable isotropic work-hardening 
coefficient. A procedure for evaluation of the material 
constants, based partly on trial-and-error curve-fitting, 
is also provided. Miller found reasonably good agreement 
of his theory with experimental results for type 30 
stainless steel at room temperature. 

This model contains both kinematic and isotropic 
hardening terms, an arrangement which leads, in Miller's 
view, to a rather complex set of equations. Miller and 



Sherby [115] simplified this in 1978 by eliminating 
kinematic hardening and allowing only isotropic hardening 
to be modelled. While this representation gave fairly 
realistic simulations for pure aluminum, it became apparent 
that an additional strengthening term was needed for 
modelling of solute-strengthened materials. The solute 
strengthening term, which is not an internal state 
variable, was added to the drag stress in the inelastic 
strain rate equation. 

In 1980, Miller noted that several existing unified 
constitutive theories, Including his own, predict 
"over-square" hysteresis loops with respect to 
experimentally observed behavior [116]. To eliminate this 
discrepancy in Miller's model, the work-hardening 

i 

coefficient in the back stress growth law, previously a 
constant, Is made an exponential function of the back 
stress and the direction of the effective stress. The 
predicted hysteresis loops using the improved equations 
compare more favorably with experimental results. 

Schmidt and Miller further improved the model by 
introducing an additional solute strengthening term 
[117,118]. This interactive term is multiplied by the drag 
stress to accompany the already present non- interact iv e 
solute strengthening term (see reference [115]) which is 


added to this product. The non-interactive term models the 



influence of 3 o 1 u t e s on yield strength while the 
interactive term models the influence of solutes on strain 
hardening. These improvements allow the model to be used 
over a broader range of loadings where solute effects are 
significant. 

Miller and his colleagues have devoted much effort to 
research in mierostructural topics which have led to the 
particular formulation and features of the model: the 
mechanisms of solutes and their effects on the yield 
strength and strain hardening behavior of a wide variety of 
alloys [119], the physical factors governing power law 
creep as a basis for the form of the inelastic strain rate 
equation [120], and experimental cyclic torsion data for 
support of many assumptions and predictions of the model 
with respect to steady-state back stresses and cyclic 
strain hardening/softening [121]. 

Recent work includes extending the theory to multiaxial 
form by Miller [122], and in a different manner by Kagawa 
and Asada [123]. In addition, It has been reported [91] 
that the latest form of Miller's theory contains four 
internal state variables in order to model the long range 
and short range components of the drag stress and the back 
stress. This form of the model may be useful in modelling 
transient temperature response [91]. 
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General Theory 

In attempting to model a broad spectrum of deformation 
phenomena, Miller uses a combination of micromechanics and 
creep phenomenology to derive the constitutive equations of 
his model [90]. It is postulated that all Inelastic 
deformation can be characterized adequately by 

e 1 = f[(o - R)/DJ, (97) 

. I 

where e is the inelastic strain rate, f is the inelastic 
strain rate function, R is the "rest stress" (or back 
stress), and D is the drag stress. To complete the model, 
it remains to determine the best possible specific function 
for equation (97) and the evolution equations for R and D. 

In stressing the need to model steady-state creep rates 
accurately, Miller uses the Garofalo hyperbolic sine 
relation for steady-state creep as a basis for deriving the 
inelastic strain rate function: 

e J = B’ [sinh(Ao) ]°, (98) 

5 3 o S 

where the subscripts "ss" refer to the steady-state 
condition, B' is a temperature dependent parameter, and A 
and n are temperature-independent constants. 

Since equation (97) must reduce to (98) for the special 
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case of steady-state creep, there must exist a function f 


1 


which causes the argument (o - R )/D to reduce to o , 

ss ssss S3 

i . e . 


f,[(o - R )/D ]=o 

1 S3 S3 SS S3 


(99) 


S3 


When this is true, equation (99) can be substituted for o 
in equation (98) to give e 1 as a function not only of 

S 3 

o , but also of R and D . Since steady-state creep is 

3 3 3 3 3 S 

only a special case, the subscripts "ss" can be dropped to 
give 


e 1 - B' (sinh{Af 1 [ ( o - R)/D]} ) n . 


( 100 ) 


By examination of warm-working data for several metals, the 
function f 1 is determined to be 


f 1 = \ C(o - RJ/D] 1 ' 5 


( 101 ) 


Substituting this into equation (100) will result in the 
explicit equation governing the inelastic strain rate: 


e 1 = B' ( sinh { [ ( o - R)/D] 1 * 5 } ) n . 


( 102 ) 


The entire temperature dependence of the model lies in 
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the constant B', which Is related to the activation energy 

for plastic flow. Below 0.6T , where T is the absolute 

m m 

melting temperature, the activation energy decreases 
linearly with temperature and B • is given as 

B' = B exp{ [-Q/0 . 6kT ] [ln(0.6T /T) + 1]}, (103) 

ra m 

where B is a material constant, k is the gas constant, Q is 
the activation energy for plastic flow, and T is the 
absolute temperature. Above 0.6T m » the activation energy 
is approximately constant and B* is given as 

B' = B exp(-Q/kT). (104) 

In either case, the expression for B’ can be abbreviated as 

B ' = B 0 ' , (105) 

where 8' is the desired temperature-dependent factor. 

The evolution equations for R and D are given in the 
standard work-hardening/recovery format as 

R = H 1 t 1 - f 2 (R,T) (106) 


and 
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D - H^e 1 - f (D,T) , (107) 

where and are constants, and and are recovery 
functions. 

In order to determine the recovery functions f^ and f^. 

Miller applies equations (106) and (107) to the special case 

of steady-state creep (where R=D=0) and substitutes 

equation (98) for e 1 to obtain 

s s 

f 2 (B 3 S - T) - “,Ss ■ H 1 B ' t = lnl'(Ao ss ): n (108) 

and 


f 3 (D 33- T) - »2 ; 33 ■ H 2 B ' tsinh(Ao gs )] n . (109) 


In a manner similar to the derivation of the inelastic 
strain rate expression, equations (97) to (100), the 
evolution equations become 


R = He 1 - H B 8' [ s i nh ( A 1 R ) ] n 


( 110 ) 


and 



HJ, B 8' [sinh( A 2 D 3 ) ] n , 


D 


( 111 ) 



where and A 2 are material constants. 

For the purposes of this research, a simplified version 
of this model will be used. Miller and Sherby [115] 
simplified the model by eliminating all kinematic hardening 
terms, allowing only isotropic hardening to be modelled, 
and incorporating a solute strengthening term: 

e 1 = B e* ( sinh { [ ( o/E ) / ( D + F o . ) 1 / 2 ] 1 * 5 } ) n , (112) 

sol 

where the modulus-compensated stress o/E has replaced the 
stress o from the earlier equations for mathematical 
convenience, F gol is the solute strengthening parameter, 
and D and F sol are placed under the square-root sign to 
produce parabolic hardening, since parabolic hardening is a 
better approximation for many materials than the linear 
hardening of the earlier equations. The evolution equation 
for D becomes 

6 = H j e 1 j - H B 0' [sinh(A 3 D 1 * 5 ) ] n , (113) 

where H and A are material constants. 

It should be noted that a slight modification has been 
introduced by this author to equation (113)* In order to 
be able to simulate cyclic loading, the inelastic strain 
rate term in equation ( 113 ) has been placed in absolute 
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value signs. This follows the practice Introduced by 
Miller, e t a 1 . , in later works [90,116-118]. 

The variable F ^ is introduced into the equations in 
order to characterize the effects of solute strengthening 
on material behavior that are found in many solute 
strengthened alloys. These effects Include plateaus in the 
yield strength versus temperature curve, local maxima and 
minima in the strain rate sensitivity versus temperature 
curve, and peaks in the apparent activation energy. By 
judicious selection of the functional dependence of F gol 
upon the other variables, one can ensure that the resulting 
equations simulate these effects. 

At low temperatures, solute atoms are so immobile that 
they cannot keep up with the moving dislocations and hence 
cannot exert any particularly strong drag force. At high 
temperatures, the solute atoms are so mobile that they can 
move easily with the dislocations and are again unable to 
exert a strong drag force. At intermediate temperatures, 
however, the average velocities of solute atoms and 
dislocations are approximately equal, and the two can 
interact strongly, producing a maximum strengthening 
effect. From these physical considerations, one can deduce 
that F go ^ should have a low value at low temperatures, 
should increase with temperature and pass through a 
maximum, and should fall to a low value at high 
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temperatures . 

Furthermore, experimental evidence indicates that 
strain rate sensitivity occurs in a fashion similar to 
temperature sensitivity. 

The particular statistical function which satisfies the 
mathematical requirements is 


sol 


F sol ,raax «P(-lClog 10 (Z) - log, 0 (Z Balt )]/2) ) 


* HlO)' 7 exp(-{[log ]0 (Z) - 303/10) 2 ), OH) 


l • I l 

where Z is e /6', F . and Z are constants, and 

the second term on the right hand side of the equation is 

included to ensure a non-zero value of F , for conditions 

sol 

when Z>>>Z 

max 

Finally, the material constants to be determined are A, 

B, n, H, F . , and Z . In addition, T and Q must be 

sol , max max m 

determined for the temperature-dependent factor 0’. 


Evaluation of Material Constants 


Q and T 

m 

The values of the absolute melting temperature T m and 
the activation energy for plastic flow for the material of 
interest may be determined from various sources in the 
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literature. With these values and by use of equation (103) 
or (10 1 !), the temperature-dependent factor 0' may be 
calculated for a given absolute temperature T. 

A 

In order to derive the relation between the constants 
A, B, and n and steady-state creep rates, one returns to 
the observation that, for the special case of steady-state 
creep, the inelastic strain rate equation (112) reduces to 
the Garofalo equation (see equation (98)): 

e* - Be' [ sinh( A o /E ) ] n , (115) 

a a sa 

• X 

where e is the steady-state creep rate and o /E has 
ss as 

replaced o . Equation (115) may be rewritten as 

3 3 

l08 io ( S3 /e,) ' 1o 8,o ( b) * 

nlog 1 0 [sinh(Ao g3 /E)]. (116) 

If this equation is an adequate representation of 
experimental data, then the constant A can be chosen such 
that log. ft ( e? 0 /6' ) is linear in log n [ sinh ( A o„ /E ) ] , as 

lUSS I U 33 


shown in Figure 1 5 . 






B and n 
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Since the value of A is now known, one can easily 
determine the constants B and n from equation (116), either 
by graphical means or by a least-squares analysis. 

F , and Z 

sol .max max 

The maximum value for F gQl is assigned arbitrarily and 

is based on an intuitive analysis of the problem of 

interest. Z is set by Miller as Z - 1(10) 10 for his 
max max 

choice of function for F 

sol 

H 

Finally, the value of the work-hardening coefficient H 
is determined on a trial-and-error basis by a best fit of 
the model through several experimental constant strain-rate 
stress-strain curves. 

In summary, Miller's theory requires some information 

about the material of interest (for T and Q), several 

m 

creep tests (constants A, B, and n), and several constant 
strain rate tensile tests (for e and H). This concludes 
the evaluation of material constants for Miller's model. 
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DETERMINATION OF MATERIAL CONSTANTS 


Each theory la evaluated and compared to experiment by 
its authors for a particular material under particular 
temperature conditions. In general, the choice of material 
and conditions varies from theory to theory. One of the 
objectives of the current research is to create a common 
experimental data base for one material at one temperature 
in order to characterize the mathematical representations 
and to compare to experiment the numerical predictions of 
each theory. 

Many of the phenomena associated with inelastic 
deformation such as rate sensitivity, the Bauschinger 
effect, and creep become significant in crystalline 
structures at elevated temperatures, that is, above three 
tenths of the absolute melting temperature. At the start 
of this research it was intended to perform testing on the 
nickel-based superalloy IN718 at 1350°F. Testing was 
restricted to room temperature conditions, however, since 
extensive difficulties were encountered with the available 
equipment with testing at high temperatures. This 
restriction led to the choice of the aluminum alloy 5086 as 
the candidate material. This alloy displays many of the 
phenomena of inelastic deformation and is already above 
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three tenths of its melting temperature at room 
temperature. Its composition and some properties are 
1 isted in Table 1 . 

A review of the material test requirements of each 
theory is shown in Table 2. It should be noted that the 
stress-drop tests are performed during steady-state creep 
and, as such, can be combined with the creep tests. The 
constant strain rate tensile tests were performed on an MTS 
810 hydraulic material testing system; the creep tests and 
the stress-drop tests were performed on a constant- load 
creep frame which was designed and constructed at Texas A&M 
University. Data acquisition in both cases was by 
graphical means, although a digital system was also 
available for the creep frame. It was decided to use a 
constant- load creep frame rather than a constant-stress 
creep frame due to difficulties associated with accurately 
determining the applied stress when using the latter 
system. As can be seen by the creep response below, the 
effects of this choice on the results are negligible. 

The results of the constant strain rate tensile tests 
can be seen in Figure 16 for four applied strain rates 

- 7 _ i -1 

ranging from 4(10) sec to 4 ( 1 0 ) sec . It can be 
seen that the material exhibits some strain rate 


sensitivity and considerable work-hardening. 
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Table 1 

Composition and Material Properties of A1 5086 


Material: A1 5086 Hill 

Composition: 0.15$Cr, 0.110$ Mn , 4.0$ Mg, remainder A1 

Melting temperature : 858K 

Elastic modulus E: 10.313(10)^ ksi 

Ultimate stress o : 40 ksi 

u 

Activation energy for plastic flow Q (at room temp.): 

27,500 cal /mole 

Form of specimens: uniaxial bar with ^-inch diameter 

gauge section and threaded ends. 
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Table 2 

Review of Material Test Requirements 


Type of test 

Kr 1 eg 

Bodner 

Miller 

Total 

Constant strain rate 
tensile test 

1 

3-4 

3-4 

3-4 

Creep 

1 

3-4 

3-4 

3-4 

Stress-drop 

3-^ 

- 

- 

3-4 



S'2 



Strain (in/in) *10 

Fig. 16 Experimental constant strain rate tensile test 


response 
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Tne results of the constant- load creep tests can be 
seen In Figure 17 for five applied stresses ranging from 
32.725 ksi to 38.850 ksi. Pertinent results of the 
constant strain rate tensile tests and the creep tests, as 
well as the stress-drop tests, are summarized in Table 3. 

Material Constants for Krieg, et al . 

In determining the material constants for this and the 
other theories under consideration, the procedures 
discussed in a previous section are followed. The reader 
may consult that section for greater detail. 

The values of c 1 and c 2 are found by applying a linear 

curve-fit to a graph of InCe 1 ) v. lnjo - a| data, where 

. 1 

e is the strain rate immediately preceding the 
stress-drop, 0 is the applied stress, and a is the back 
stress. This is shown in Figure 18; the slope of the 
linear fit is c 2 and the vertical intercept is ln(c 1 ). In 
this manner, c 1 is calculated to be 7.177 and c 2 is 
1 . 459(1 0)" 9 . 

Under steady-state creep, the growth law for the back 
stress becomes (equation (71) repeated) 

(c^/c^) = [ a 2 / ( c 1 j q j C 2 ) ] [exp(c 5 a 2 ) ~ 1]. (71) 


Applying this equation to two experimental cases results in 
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Table 3 

Pertinent Results of Experimental Tests 

Constant strain rate tensile tests 
Elastic modulus = 10.313(10 ) 3 ksi 
Elastic limit = 28.000 ksi 


0.2$ offset yield stress = 30.250 

ksi 

for ^(lO)" 4 

-i 

sec 

30.000 

ksi 

4(1 0)' 5 

-1 

sec 

29.750 

ksi 

4(1 0)" 6 

-1 

sec 

29.500 

ksi 

1 

o 

w 

*=r 

-1 

sec 


Creep tests 

Steady-state creep rate *» 


1 .205(1 0)' 8 

sec 1 

for 32.725 

ksi 

4 . 1 67 ( 1 0 ) ”" 8 

sec 1 

35.433 

ksi 

3.1 25(1 0)' 8 

-1 

sec 

36. 437 

ksi 

4. 435 (1 0 )' 8 

sec 1 

36.705 

ksi 

4.546(1 0)' 8 

-1 

sec 

38.850 

ksi 


Stress- drop tests 

Strain rate before applied stress back stress 

stress drop 


0 . 1 000 ( 1 0 ) " 7 

-1 

sec 

3 4 . 402 

ksi 

33-107 

ksi 

0.4243(1 0)" 7 

-1 

sec 

35.433 

ksi 

33-804 

ksi 

2.658(1 0)" 7 

-1 

sec 

36.705 

ksi 

34.659 

ksi 





two nonlinear equations in two unknowns which can be solved 
for (c^/Cjj) and by iteration. Solving for the cases 
with the highest and lowest values of a (see Table 3) 
yields ( c _ / c ) = 3.255(10) 25 ksi 2 sec and c c = 0.03030 
ksi -2 . 

To determine , a point is chosen on the primary creep 
curve where the creep rate is still twice the eventual 
steady-state rate. The effective stress £ at this point is 
calculated by inverting the Inelastic strain rate equation 
(68) (equation (73) repeated here): 

C = (e I /c 1 ) WC *. (73) 

The value of the back stress is then given by a = o - £. 
Using the creep test with o = 35.433 ksi gives the 
following information: 

o i = 35.433 ksi, 

ej = 8.3330 0) -8 sec” 1 , 

t 1 *= 450 sec. 

This gives ^ = 1.717 ksi and = 33.716 ksi, where the 
subscript 1 refers to the point on the primary creep curve. 
The growth law for a, equation (69), can be integrated 
numerically by separating variables (equation (72) 


repeated ) : 
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Vi - C VK - “!° 2 - 

(c^/c^a 2 [exp(c^a^) - 1]} ^a. (72) 

g 

With a 1 and o 1 known, c^t is calculated as 2.126(10) ksi 
sec. From this, c^ = 5.185(10) ksi and 
c ^ = 1 .593(1 0)~ 21 ksi -1 sec" 1 . 

This concludes the material parameter determination for 
the theory of Krieg, e t al . 


Material Constants for Bodner, et al . 

For the present strain rate range of interest, D Q is 
i| - 1 

set as 1(10) sec 

To determine n and , a saturation stress must be 

determined analytically since the experimental 

stress-strain curves are not saturated (see Figure 16). 

This is accomplished by extrapolating to the point of zero 

P 

Y on the Y v. o graph, where Y=do/dW and is given by 

P P 

differentiating a curve-fit of W = W (o). It is found 

P P 

that quadratics provide the best curve-fit for W = W (o): 
W P = 4.679 - 0 . 3 679 o + 0 . 7 46 2 ( 1 0 ) " 2 o 2 (117) 

for e = 4(10) ^ sec 1 , and 
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Values of c 1 are calculated from e 1 = e (1 - 

E Ac 

and the resulting values of Z are shown in Table 4 below. 


Table 4 . Calculated Values of Z used 


• , - 1 \ 
e (sec . ) 

in the 
o (ksi) 

Determination of 
t 1 (sec } ) 

Z 0 

Z (ksi) 

r>- 

1 

o 

«=r 

28.250 

3. 6474(1 0 )~ 7 

34.102 

4 ( 1 0 ) ~ 6 

28.375 

3.5152(10)“ 6 

34.090 

-fcr 

o 

1 

V7I 

28.500 

3.3536(1 0)~ 5 

34.060 

4(1 O)" 4 

28.688 

3.2243(1 0)"** 

34.081 


An average value of 34.00 ksi is taken for Z Q . 

The value of , the minimum value of Z, is taken to be 
the lowest value of Z for steady-state creep, although one 
may simply set Z 2 « Z Q . The latter is chosen, giving 
Z 2 = 33.50 ksi. 

The values of the recovery coefficient A and the 
recovery exponent r are determined by calculating a value 
of Z for the steady-state region of each creep test by 
using equation (94) (see above). Calculated values of Z 
are shown in Table 5 below. 



Table 5. Calculated Values of Z U3ed 


in the Determination of A and r 


°ss 
( ksi ) 

• I 
e ss 

(sec 1 ) 

Z 

( ksi ) 

ln[m (Z 


ln[(Z-Z 2 )/Z 1 ] 

32.725 

1 .205(1 0)“ 8 

39.758 

-2.091 


-12.38 

36.437 

3.1 25(1 0)" 8 

44.1 92 

-1 .520 


-12.37 

36.705 

4. 435(1 0)" 8 

44.489 

-1 .491 


-12.14 

These values of Z are 

then 

used to 

construct a 

ln[m ( Z 1 

- Z) W P ] v. 

ln[ ( Z 

- z 2 )/z, 

] graph, 

as shown in 

Figure 21 

. The slope 

of a 

linear curve-fit 

for this graph 


is r; the vertical intercept is ln(AZ^). In this manner, A 
is calculated to be 1.447(10) ^sec 1 and r is 0.2290. 

This concludes the material constant determination for 
the theory of Bodner, et al . 
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® Data points 

Best fit 
line 


ln[m(Z^-Z)W P ] 



Fig. 21 Graph of in[m (Z^Z) W P ] v. ln[(Z-Z 2 )/Z ] 
for Bodner’s Model 
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Material Constants for Miller, e t al . 

The values of the activation energy for plastic flow Q 

and the melting temperature T^ can be found in a materials 

handbook. For Al 5086 at room temperature, Q is 27500 

cal/mole and T is 858K. 
m 

Since room temperature is less than 0.6T for Al 5086, 

m 

equation (103), repeated below, is used to calculate the 
temperature-dependent factor e*: ' 

0* = exp { [-Q/ 0 . 6kT ] [ In ( 0 . 6T /T) + 1]), (103) 

ra m 

where T and are absolute temperatures. Thi3 gives 
0 * as 8 . 7.4(1 0)' 1 9 . 

The value of A is chosen such that a graph of 

log, .( e 1 /0’) v. log, A [sinh( Ao /E)] for steady-state 
loss in ss 

creep data is linear. This is shown below in Table 6 and 
in graphical form Figure 22 for several values of A. 



Table 6. Calculations for Various Values of A 


1Og 1O U S3 /0,) °ss l°8 10 [sinh(Ao 33 /E)] 


(sec 1 ) 


( ksi ) 

A = 50 A = 200 A 

= 1000 

1 .205(1 0)" 8 

10.14 

32.725 

-0.7977 -0.1687 

1 .076 

4.1 67(1 0)” 8 

10.68 

35.443 

-0.7628 -0.1292 

1.191 

3.1 25(1 0)" 8 

10.55 

36. 437 

-0.7506 -0.1 1 53 

1 .233 

4. 435(1 0)” 8 

10.71 

36.705 

-0.7474 -0.1116 

1 .244 

4.546(1 0)” 8 

10.72 

38.850 

-0.7225 -0.0826 

1 .335 


As can be seen in Figure 22, t ^ curve’ is near-iy linear for 
A = 50, with no significant improvement for A < 50. 

With this value of A, the values of B and n can be 
determined by constructing a linear curve-fit for the data 
in Figure 22. The slope of the curve-fit is n; the 
vertical intercept is log 1Q (B). In this manner, n is 
calculated to be 7.701 and B is 2.410O0) 18 sec . 

The maximum value of F gol is set arbitrarily as 

- 7 

F , = 1(10) . Simulations of experimental tests have 

sol , max 

shown that the choice of this value does not influence the 
results significantly. 

The work-hardening coefficient H is set arbitrarily to 
give the best fit of simulations of constant strain rate 
tensile tests to experimental data. This value is set as 


7.5(10)”** ksi. 




Finally, the initial value of the drag stress, 

is calculated by inverting the inelastic strain rate 

equation . (1 1 2) and using values of e* and <? at yield. In 

- 1 ; 

this manner, D^ ni is calculated to be 3-386(10) 

This concludes the material constant determination for 
the theory of Miller, et al . 

A sun ,ary of the material constants for each theorv is 


given in Table 7 



Table 7. Summary of Material Constanta for all Models 


Krieg, etal ; E = 10.313(10) 3 ksi 



C 1 

= 7.177 



C 2 

= 1 . 459 (1 0)" 9 



C 3 

= 5.1 85(1 0) ^ 

ksi 


c 4 

= 1 .593(1 0)" 21 ksi -1 sec' 1 


°5 

- 0.03030 ksi 

-2 

Bodner, et al. 

E - 

10 . 313 ( 1 o) 3 

ksi 


n = 

10.30 



Z 1 

= 46.59 ksi 



m = 

1 . 558 ksi -1 



z o 

= 34.00 ksi 



Z 2 

= 33.50 ksi 



A = 

1 . 447 (1 0 )“ 7 

-1 

sec 


r = 

0.2290 


Miller , et al . 

E = 

1 0 . 3 i 3 ( 1 0 ) 3 

ksi 


e ' 

- 8.74(1 0)' 1 9 



A = 

50 



n = 

7.701 



/ 

B = 

2.41 0(1 0) 1 6 

-1 

sec 


H = 7 . 5(1 0)~ H ksi 

D i n i t = 3.386(10) 



COMPARISON OF THEORY TO EXPERIMENT 


The third and final object tve of this research is to 
compare the numerical predictions of each theory to 
experiment by integrating the constitutive equations 
subject to specific load histories. Before proceeding to 
the quantitative comparisons, however, there follows a 
brief qualitative review of the predictive capabilities of 
each theory. 

Qualitative Review of Predictive Capabilities 
All the theories under consideration are similar in 
that they are unified theories; that is, the rate-dependent 
creep and the rate- independent plasticity components have 
been combined into one inelastic strain term. Each theory, 
in the form under consideration, is isothermal in that it 
is not able to model transient temperature response. 
Furthermore, each model is able to exhibit strain rate 
sensitivity through the exponent in the inelastic strain 
rate equation. 

The theory of Krieg, et al . , contains two internal 
state variables: the inelastic strain and the back stress. 
The back stress represents kinematic hardening, which 
allows the Bauschinger effect to be modelled. Krieg, 



e t a 1 . , have assumed the hardening function in the back 
stress growth law to be constant, which may cause 
"over-square” hysteresis loops. This may be alleviated by 
introducing a dependence of the hardening function on the 
current value of the stress or back stress. Due to the 
lack of a representat ion for isotropic hardening (the drag 
stress), this theory may not be able to model cyclic strain 
hardening. Finally, it can be seen that anelasticity 
(strain recovery at zero load) can be predicted by this 
theory by setting o = 0 in the inelastic strain rate 
equation (68); relaxation is predicted by setting e =0; 
and creep is predicted by setting o ■= 0 . 

The theory of Bodner, et al . , also contains two 
internal state variables: the inelastic strain and a terra 
representing isotropic hardening. The isotropic hardening 
growth law contains dynamic as well as thermal recovery. 

The presence of isotropic hardening allows cyclic strain 
hardening to be modelled, although the absence of kinematic 
hardening may preclude the ability to model the Bauschinger 
effect. Finally, it can be seen that anelasticity cannot 

. i 

be modelled since e is zero when setting o = 0 in the 
inelastic strain rate equation (86); creep and relaxation 
can be predicted, however, by setting o = 0 and e = 0, 
respect ivel y . 


The theory of Miller, et al., 


contains two internal 
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state variables: the Inelastic strain and the drag stress. 
The drag stress represents isotropic hardening, which 
allows cyclic strain hardening to be modelled. This model 
is similar to that of Krieg, et al . , in that the hardening 
function in the drag stress growth law is assumed to be 
constant, which produces "over-square" hysteresis loops. 

The model is similar to that of Bodner, e t al . , in that the 
lack of a back stress may preclude the ability to model the 
Bauschinger effect as well as anelas t i c i ty . Finally, creep 
and relaxation are modelled in a manner similar to that of 
the other two theories. 

Review of Integration Technique 
Each theory is represented mathematically by a set of 
differential equations. Since these differential equations 
are numerically "stiff", causing higher-order integration 
schemes to become unstable [124], a first-order forward 
integration scheme with small, variable time steps was 
used. For example, Krieg's Inelastic strain rate equation 
((68) repeated here) 

e = c^ j o - a | C 2 sgn(o - a) (68) 


becomes 



£ 1 ( t + At) = c 1 |o(t) - a(t) | C 2 A t + e 1 ( t ) (119) 

where t represents an initial state and t + At represents 
the incremented state. Similarly, the growth law for the 

back stress (equation (69)) becomes 

X 

a ( t + At) - {c 3 [ ( e 1 ( t + At) - e I (t))/At] - 

c^a 2 ( t ) [ exp (c 5 a 2 ( t ) ) - 1 ] sgn ( a ( t ) ) } A t + a(t).(120) 

Finally, the stress is updated by 

o ( t + At) = E { e - [ ( e 1 ( t + At) - e 1 ( t ) ) / At ] } At + 

O(t). (121) 

This concludes the review of the integration technique used 
in this research . 

Comparison of Numerical Predictions to Experiment 

In this section, the ability of each model to predict 

the response to a specific load history is analyzed. The 

first experimental test to be modelled is a constant strain 

• 5 * * 1 

rate tensile test with e = 4(10) sec , shown in 


Figure 23- 







Krieg, e t a 1 . : Constant Strain Rate Tensile Test 


It can be seen from Figure 24 that this theory does not 
compare favorably with experiment for this constant strain 
rate tensile test. It Is believed by this author that this 
Is due largely to the method of determining the material 
constants. As discussed in an earlier section, the ability 
to determine a zero creep rate or merely the absence of a 
resolvable creep rate in a stress-drop test depends greatly 
on the resolution of the data acquisition equipment and is 
extremely difficult and time consuming at best. A 
parameter variation study showed slight improvements in the 
stress-strain curve in Figure 24 but led to computational 
difficulties in simulating cyclic load histories and 
complex load histories. It can be concluded from this that 
improvements in the measurement of back stresses are 
required. For example, the methods of Blum and Finkel 
[98], in which the instantaneous strain change and the 
maximum strain change after unloading are compared, and of 
Walker [58], in which one applies hold times on the 
unloading side of a saturated hysteresis loop and 
extrapolates to the point of zero relaxation, merit further 
study. 

Bodner, et al . : Constant Strain Rate Tensile Test 

The theory of Bodner, etal . , compares favorably to 



Stress (KSI) 

2.50 7.50 12.50 17.50 22.50 27.50 32.50 37.50 42.50 47.50 




experiment, as can be seen in Figure 25. Work-hardening is 
predicted well although the el ast ic- inelast ic behavior 
transition is too sharp. It is found that this is due to 
the theory predicting values of Z that are too high at low 
stresses, possibly requiring a modification of the 
hardening term in the growth law for Z. 

Miller, et al . ; Constant Strain Rate Tensile Test 

In determining the material constants for this theory, 
the strain rate exponent n was calculated to be 7.7, 
resulting in the response shown in Figure 26. A parameter 
variation study indicated that a value of n « 6 appeared to 
fit the experimental data much better. This indicates a 
problem inherent in many constitutive theories: the 
determination of material constants is not rigorous and the 
difficulties are compounded by the fact that many constants 
are interdependent and, in fact, may be assigned arbitrary 
values. In the interests of further numerical comparisons, 
a value of n = 6 is chosen for this model. With this 
choice, this theory predicts the yield point and 
work- hardening fairly well. 

Finally, a comparison of all the models for a constant 
strain rate tensile test can be seen in Figure 27- 


a EXP» STRAIN RATE » 4E-5 
O 0OD>CR. STRAIN RATE - 4£-5 



Fig. 25 Constant strain rate tensile test. Comparison of 
theory of Bodner, et al . , to experiment. 
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Fig. 26 Constant strain rate tensile test. Comparison of 


theory of Miller 


e t al 


• I 


to experiment 
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Fig. 27 Constant strain rate tensile test. Comparison of 


all theories to experiment. 



Many theories use creep tests or constant strain rate 
tensile tests as information for the material parameter 
determination; consequently, one would expect fairly 
accurate simulations of these types of load histories. A 
more rigorous comparison, however is a cyclic load history, 

• — ij - 1 

as shown in Figure 28 with e = _+ 1 ( 1 0 ) sec for cycles 
1-10. The experimental data exhibit a pronounced 
Bauschinger effect and considerable cyclic strain 
hardening, although it appears that the hysteresis loop 
saturates fairly quickly. 

Krieg, et al ♦ : Cyclic Loading Test 

The cyclic response of the theory of Krieg, et al . , is 
shown in Figure 29 for cycle 1 . It can be seen that the 
theory does indeed predict a Bauschinger effect, although 
it is not as pronounced as the experimental data. The 
model predicts a harder material than is the actual case, 
resulting in the very thin and, as expected, "over-square" 
hysteresis loop. This may be due to the ratio c^/c^ 
(calculated from stress-drop data) being too large. 

A comparison of the prediction for cycles 1 and 10 can 
be seen in Figure 30, where, as expected, the model does 
not display significant cyclic strain hardening. However, 
since the material exhibits cyclic strain hardening, the 
difference between the saturated experimental and 
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B KRIEG. STRAIN RATE - V-1E-4. CYCLE l 
© EXPi STRAIN RATE = V-1E-4. CYCLE 1 



Fig. 29 Cyclic loading test: Cycle 1. Comparison of 
theory of Krieg, et al . , to experiment. 




Stress (KSI) 

-45. DO -35.00 -25.00 -15.00 -5.00 5.00 15.00 25.00 35.00 45.00 





1 22 


predicted loops diminishes somewhat, as shown in Figure 31 • 

Bodner, et al . ; Cyclic Loading Test 

Figure 32 shows that this theory models cycle 1 fairly 
well, although the ela s t i c- i nel a s t ic behavior transition is 
too sharp (as noted in earlier tensile data also). The 
absence of the Bauschinger effect suggests that a 
directional harden! ng term is necessary. This can also be 
seen when examining the Y v. o graph used earlier in 
calculating some of the material constants. Although the 
theory assumes a linear relationship between Y and o, the 
experimental data suggest a bilinear relationship. This is 
examined in greater detail by Bodner, et al . , in recent 
work [109-111]. The upper slope in this bilinear 
relationship governs directional hardening and the lower 
slope governs isotropic hardening. It can be seen easily 
that including such a directional hardening term could 
provide the Bauschinger effect and also decrease the height 
of the hysteresis loop. 

The aforementioned observations can also be seen in 
Figures 33 and 3^, which compare prediction for cycle 1 to 
prediction for cycle 10 and prediction to experiment for 
cycle 1.0, respectively. 
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a KRIEG. STRAIN RATE - V-1E-4. CYCLE 10 
© EXPi STRAIN RATE = -/-IE-4. CYCLE 10 



Fig. 31 Cyclic loading test: cycle 10. Comparison of 
theory of Krieg, et al . , to experiment. 




Strain <in/in) 


* 10 '' 


Fig. 32 Cyclic loading test: cycle 1 . Comparison of 
theory of Bodner, elal. , to experiment. 









Stress (KSI) 
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0 BODJCfi. STRAIN RATE - -/-IE-4, CYCLE 10 
© EXP. STRAIN RATE - V-1E-4. CYCLE 10 



Fig. 3^ Cyclic loading Lest: cycle 10. Comparison of 


theory of Bodner, et al. , 


to experiment. 




Miller, e t a 1 . : Cyclic Loading Test 

The similarities between this model and that of Bodner, 
e t al . , can be seen in Figure 35 in that the theory does 
not model the Bauschinger effect due to the absence of a 
directional hardening term. It can also be seen that the 
assumption of a constant hardening function for the drag 
stress causes an "over-square" loop, although this is not 
as pronounced here as with the other two models. 

Comparing the predictions of cycle 10 and cycle 1 shows 
extensive cyclic strain hardening (see Figure 36). This is 
much more apparent in the theory- to-exper iment comparison 
for cycle 10, as shown in Figure 37, and suggests a need 
for a revision of the hardening function. 

Finally, a comparison of all models for cycles 1 and 10 
is shown for completion in Figures 38 and 39, respectively. 
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Fig. 35 Cyclic loading test: cycle 1 . Comparison of 


theory of Miller, e t al . , 


to experiment. 
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Fig. 36 Cyclic loading test: cycles 1 and 10. Prediction 


or theory of Miller, et al . 
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Fig. 38 Cyclic loading test: cycle 1 . Comparison of all 


theories to experiment 
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Fig. 39 Cyclic loading test: cycle 10. Comparison of all 


theories to experiment. 



The final set of comparisons involves the prediction of 
the response to a complex loading history. Many models are 
designed specifically for creep response or cyclic response 
with somewhat limited capability for a general, complex 
response. The loading history is shown in Figure 40; the 
experimental response is shown in Figure 41. 

Krieg, et al . : Complex Loading History 

Figure 42 shows that this model over-pred lets the 
stress, indicating that the Inelastic strain is growing too 
slowly. This implies that the back stress is growing too 
quickly and that the hardening/recovery ratio is too large. 
Further evidence of this is demonstrated by the 
considerable amount of relaxation in compression. As 
before, this discrepancy may be due to the method of 
material constant determination. 

Bodner, et al . : Complex Loading History 

Figure 43 shows that this theory predicts the complex 
loading response very well. The model appears to predict a 
harder material with less relaxation than the experiment, 
again suggesting that the hardness Z is too large. A large 
Z causes a large e, leading to an over-prediction of 
stress. Similarly, during relaxation (when e = 0), a 
large Z causes a small o, leading to little relaxation. 



Fig. 14 0 Experimental complex loading history 
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Fig. iJ2 Complex loading test. Comparison of theory of 
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Fig. 43 Complex loading test. Comparison of theory of 
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Miller, e t a 1 . : Complex Loading History 


Figure 44 shows that this model o v er- pr ed i c t s the stress 

during loading and unloading. This may be due to the fact 

that the hardening constant was set arbitrarily by fitting 

theory to experiment for a given strain rate, whereas the 

complex loading test involves strain rates upto two orders 

of magnitude higher. The large amount of relaxation is 

then brought on by the large value of stress. 

Furthermore, it should be noted that, at approximately 

40 seconds, the model experiences a negative stress rate 

immediately after the applied total strain rate is reduced. 

Again, this is due to the large value of stress causing the 

inelastic strain rate momentarily to be larger than the 

. T 

total strain rate. This decreases o until e is less than 
e, at which point the response returns to normal. 

A solution to this could be an improved method for 
determining H, as well as a more complex hardening law for 
the drag stress. 

Finally, a comparison of all theories to experiment for 
the complex loading history is shown in Figure 45 for 


completion. 
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SUMMARY AND CONCLUSIONS 




The objective of the research in this thesis has been 
to review and implement three theories for the prediction 
of inelastic deformation in crystalline structures. These 
theories are based on considerations of the raicrostructural 
behavior of the material and are the theories of Krieg, 
e t a 1 . , Bodner, et al . , and Miller, et al . 

The thesis opens with a review of the mechanics of 
continua with internal state variables and a review of the 
historical development of constitutive modelling. It was 
shown that the concept of internal state variables is a 
suitable framework for comparison of all the models 
discussed herein and that the primary differences in these 
models are only in the number of internal state variables 
chosen and the growth laws used to obtain these variables. 

A detailed discussion of each model chosen for 
comparison follows, in which the chronological development 
of the model is given and the theory on which the model is 
based is reviewed. This section closes with a description 
of the method of determination of material parameters and 
an assessment of the material test requirements of each 
model . 


In order to perform a numerical comparison of the 



models, the aluminum alloy A1 5086 was chosen as the 
candidate material. An experimental data base was 
established, from which the material constants for each 
model were calculated. 

Finally, the uniaxial mathematical representations of 
each model have been encoded and implemented in computer 
programs in order to simulate various load histories. 

These simulations are compared to experimental data and 
comments are made on their accuracy. 

In conclusion, it has been demonstrated that the 
theories discussed herein are suitable models for the 
prediction of inelastic deformation in crystalline 
solids. It Is also apparent, however, that further 
research is necessary in several areas in the development 
of these theories. 

For example, the theories lack, to varying degrees, 
clear, concise, and rigorous methods to determine material 
constants from simple experiments. The theory of Krleg, 
etal . , uses the somewhat controversial and difficult to 
implement stress-drop test, although the calculation of 
constants is simple. This dependence on stress-drop test 
results may adversely affect the performance of this model. 
The theory of Miller, e t al . , uses standard creep tests, 
but values of several constants are set arbitrarily to 


obtain an optimum curve-fit. Although this is also true to 



a lesser degree for the theory of Bodner, et al . , it is 
found that this theory provides the best balance between 
ease of testing and determination of material constants for 
the theories discussed herein. In addition, it is found 
that all three theories exhibit substantial material 
parameter sensitivity to small variations in input. 

It can be seen in the previous section that the models 

in general compare favorably to experiment. However, it is 

clear that the full representations of the models, 

incorporating both isotropic and directional hardening 

terms, may be more appropriate for modelling phenomena such 

as the Bauschinger effect and cyclic strain hardening than 

the simplified versions examined herein. In addition, 

these theories are proposed mainly for intermediate to high 

temperature conditions and may have failed to model fully 

the phenomena at the temperature used in this research 

(0.3T ). 
m 

Another area of improvement is in the use of internal 
state variables. Many hardening functions are taken to be 
constants, whereas functions of the current stress level or 
back stress level may be warranted to prevent an 
overprediction of stress and an "over-square" transition 
from elastic to inelastic behavior. 

Based on these assessments, it is felt that the theory 
of Bodner, et al., provides the best balance between ease 



i im 


of implementation and accuracy of results for the material 
and test conditions of interest for this research. The 
theory of Krieg, et al . , is affected adversely by the 
complexity of the required material tests, while the theory 
of Miller, et al . , is affected adversely by the lack of 
rigor in calculating constants. 

As a result of this research, it is the experience of 
this researcher that data acquisition and control equipment 
of the highest resolution and accuracy is an absolute 
necessity to ensure credible experimental results, 
particularly in specialized tests such as the stress-drop 
test . 

Finally, it is hoped that the research presented in 
this thesis may answer questions relating to the modelling 
of inelastic deformation in crystalline solids and lead to 
further research in this field. 
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